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ABSTRACT 

This  thesis  describes  the  development  of  a  computer  program  that 
calculates  the  propagation  loss  for  low  frequencies  in  a  shallow  ocean 
given  the  depth  of  source  and  receiver,  the  sound  speed  profile  of  the 
water,  the  frequency  of  the  source,  and  the  impedance  and  sound  speed  in 
the  bottom.  The  program  does  this,  by  computing  the  sum  of  normal  modes 
for  a  specified  set  of  boundary  conditions.  At  the  surface,  perfect 
pressure  release  is  assumed,  and  the  boundary  condition  at  the  bottom  is 
one  of  impedance  mismatch.  An  effort  was  made  to  develop  a  Fast  Field 
Program,  which  would  use  a  FFT  to  predict  propagation  loss  at  a  variety 
of  ranges  by  solving  for  a  discrete  set  of  wave  numbers,  but  development 
was   not   completed. 
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I.    INTRODUCTION 

A  wide  variety  of  methods  are  available  for  determining  the  pressure 
field  from  a  point  source  in  a  stratified  ocean  of  constant  depth. 
Figure  1.1  shows  that,  by  using  various  transforms,  it  is  possible  to 
derive  expressions  describing  the  pressure  field  in  terms  of  the  linear 
wave  equation,  the  Helmholtz  equation  and,  finally,  the  field  integral. 
Methods  associated  with  these  expresssions  vary,  and  each  has  relative 
advantages  and  disadvantages  in  terms  of  speed,  accuracy,  and  ability  to 
deal  with  complex  situations.  Unfortunately,  no  single  method  is  supe- 
rior in  all  respects  and  there  must  be  trade-offs.  In  some  situations, 
accuracy  might  be  sacrificed  for  speed,  or  the  ability  to  deal  with 
unusual   boundaries  may  override  the  importance  of  speed  or  accuracy. 

Figure  1.2  shows  the  common  methods  of  solving  the  field  integral  to 
determine  the  pressure.  The  Normal  Mode  Solution,  the  Method  of  Stee- 
pest Descent,  The  Multiple  Scattering  Method  and  Direct  Numerical  Integ- 
ration are  "exact"  methods.  The  Fast  Field  Program  (FFP),  a  method  of 
numerical  integration,  must  be  considered  an  approximation,  not  only 
because  it  approximates  a  Hankel  function  by  an  exponential,  but  also 
because  of  its  reliance  on  a  finite  number  of  samples.  Since  Stationary 
Phase  normally  requires  approximations,  then  so  does  its  derivative,  the 
Method  of  Images.  The  Wentzel-Kramers-Brillouin  (WKB)  Integral  normally 
involves  approximations  so  both  Ray  Theory  and  the  WKB  Mode  Equation 
will    seldom    be    exact.      It    is    possible    to    find   exact    solutions    using   the 
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Normal  Mode  Technique,  but  in  this  paper,  as  will  be  seen  in  Chapter  II, 
practical  assumptions  relegate  the  solutions  to  the  approximate  cate- 
gory. This  thesis  is  limited  to  consideration  of  the  Normal  Mode  Solu- 
tion, which  is  computed  by  the  program  EXACT,  and  a  Fast  Field  Program, 
computer    program   FFP. 

Chapter  II  derives  the  field  integral,  then  develops  the  Normal  Mode 
Solution,  and  describes  the  Fast  Field  Program,  FFP.  The  basic  differ- 
ence between  the  methods  involves  the  requirement  of  the  Normal  Mode 
Solution  for  an  accurate  estimation  of  the  horizontal  wave  numbers  for 
each  of  the  modes.  The  FFP  does  not  calculate  the  modal  wave  numbers 
but  computes  the  pressure  contributions  for  a  large  number  of  pre- 
defined horizontal  wave  numbers.  An  expression  to  describe  losses  and 
phase  shift  in  terras  of  the  impedance  and  sound  speed  ratios  between  the 
bottom   and  the   water   is   also   derived  in  Chapter  II. 

Chapter  III  describes  the  programs  EXACT  and  FFP.  EXACT  solves  the 
field  integral  by  first  finding  good  approximations  of  the  mode  wave 
numbers  using  a  matrix  eigenvalue  technique  and  then  refining  them.  FFP 
calculates  the  pressure  field  at  a  discrete  number  of  horizontal  wave 
numbers  without  regard  to  modal  values.  Both  programs  use  the  vector 
sum  of  the  components  of  pressure  from  the  constituent  wave  numbers  to 
calculate  the  propagation  loss  between  the  source  and  the  receiver, 
taking  into  account  absorption  in  the  water,  losses  to  the  bottom,  and 
phase   changes   at   the  surface   and   bottom. 

Chapter  IV  describes  the  testing  procedures  and  the  results.  Within 
this  chapter  the  testing  criteria  are  defined  and  the  adaptive  nature  of 
the    development    process    is    discussed.      Test    results    provided    some    very 


instructive     lessons     which    gave     an    insight    into    how    normal     modes    can 
describe   pressure   variations  in   both  the   vertical   and  horizontal. 

Chapter    V    discusses    the    results    and    some   of    the    weaknesses    of    the 
program. 


10 


Exact  Wave 
Solutions 


Approximate 
Solutions 


Geometrical 
Acoustics 


Basic  Phyics 
Cons  of  Mass 


Linear 

I 


Lin  Wave  Eqn 
I     1  d^ 


Fourier  Transform 
(t  -v  to) 


Helmholtz  Eqn 


V^P  =  K^P 


F((o) 


Hankel  Transform 


Field  Integral 
Point  Source 

F  fuv  „(  1)  .J.  s^  ,- 
7—  -77-  "  (^r)c,dfc, 
4tt   VI       o 


Numerical 
solution 


Method  of 
Characteristics 


Numerical 
Solution 


Numerical 
Finite 
Element 
Model 


Separation  of 

Variables 

P(r,z)=R(r)Z(z) 


Mormai  Mode 
Solution 


Coupled  Modes 

P  =  E'J  (r)u  (r  ,z) 
n    n 


Normal  Mode 
Sclution 


■lethod  of 
J  Steepest  Descent 


Direct  Numerical 

Solution  FFP 


Ray 
Theory 


-^ 


tikonai 
Equation 


Ray 
Theory 


Parabolic 

Equation 

1   P(r,z)    = 

1'-,  r 
A(r,z) 

/r 

WK3 


->  Approxima Lion  ! — ->  Theory 


Ray 


Perturbation 
Solutions 


Adiabi 


->    i-louo  I 

I  Sol ut ions  I 


Figure    1.1.      Methods    of    Solving    for    the    Pressure    Field 
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II.    THEORY 

In  this  discussion  of  normal  modes,  frequency  will  be  very  low  so 
the  absorptive  losses  due  to  the  water  will  be  negligible,  except  at 
long  ranges.  However,  absorption  in  the  bottom  will  often  be  very  high, 
especially  for  the  higher  modes.  The  theory  was  used  to  derive  formulae 
for  the  lossless  case;  then  absorption  was  introduced  as  an  imaginary 
component  of  the  horiztonal  wave  number.  Although  absorption  in  the 
water  is  low  for  low  frequencies,  it  proved  expedient  to  input  abnor- 
mally high  absorption  rates  into  the  test  programs.  By  doing  this,  it 
was  possible  to  verify  that  the  complex  numbers  were  producing  consis- 
tent  values. 

The  theoretical  approach  involves  a  series  of  transforms,  some  of 
which     require     approximations     for     simplification.  These     transforms 

successively  take  the  wave  equation  from  the  time  domain  to  the  fre- 
quency domain,  from  three  dimensions  to  two,  from  range  dependent  to 
wave  number  dependent,  and  finally  to  a  function  that  describes  pressure 
as  a  function  of  the  depths  of  the  receiver  and  transmitter,  the  hori- 
zontal wave  number  and  the  horizontal  range.  Normal  mode  theory 
enables  one  to  consider  the  total  pressure  as  the  sum  of  pressures  from 
many  components.  Each  mode  is  stimulated  by  the  transmitter,  or  sensed 
by  the  receiver,  to  a  degree  that  depends  on  the  positions  of  the 
transmitter  and  receiver  in  relation  to  the  modal  depth  function  curve. 
This   curve  is   described  in  Chapter  IV. 
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In  order  to  keep  the  program  simple,  several  assumptions  must  be 
made;  some  limit  the  practicality  of  this  particular  program  and  could 
be  overcome  in  a  more  comprehensive  routine,  while  others  are  of  little 
consequence. 

The  main  assumptions,  some  of  which  will  be  discussed  in  detail 
later,   are: 

1.  source    is    cw    and   monochromatic;    it    emits    a    continuous    wave 
at   a   constant   frequency; 

2.  source   is   a   point   radiating  uniformly   in   all   directions; 

3.  medium  is   homogeneous  in  all  respects  in  the  horizontal; 
U.      bottom   and  surface   are  flat   and   parallel; 

5.  surface  is   perfectly  reflecting; 

6.  bottom    can   be   completely  described   by   its   impedance; 

7.  all    energy  transmitted  into  the    bottom   is   lost,    and 

8.  branch   line  integrals   can   be   ignored. 

A.      THE    FIELD   INTEGRAL 

An  omnidirectional  point  source  operating  at  frequency  f(t)  is  posi- 
tioned at  a  depth  z  in  a  water  mass  that  is  uniformly  homogeneous  in 
sound  speed  in.  the  horizontal  plane  and  can  be  described  in  the  vertical 
plane   by  the  sound  speed,   c(z),   at    any   depth. 

As  long  as  the  variations  are  linear,  the  pressure  is  given  by  the 
acoustic  wave   equation: 


''P   -    p(^a    3^    -    -f(t)«s-3,)  2.1 


where: 

p  -  p(3,t)  pressure  as   a  function  of  position  and  time, 
f(t)   =•   function  describing  the  source  amplitude, 
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c(z)   =  depth  dependent  sound  speed,   and 
s   =   position  in  Cartesian   coordinates. 

A    Fourier    transform    from    time    to    frequency,    'jo,    corresponds    to    con- 
sidering the   equation  for   distinct   frequencies   and   yields: 


,  .2 

V^P    +    -^    P   =    -F(aj)6(s-So)  2.2 


where: 
P   =   P(s,a3)   a   function  of   position  and  frequency, 
(jj  =   frequency  in  radians/sec,   and 
F(a))   =   transform   of   the   driving  function,   f(t). 

Because  the  sound  field  is  independent  of  the  azimuthal  angle,  this 
transform  can  be  formulated  in  the  more  manageable  cylindrical  co-ordi- 
nate system.     In  the  new   system   the  wave   equation  appears   as: 


3^P    .     1     3P         8^P        ,  ,^  r..    s    ^^^z-Zo) 


+    k^P    =    -F(oa)    — :r dr  2.3 


dr^        r    dr         3z^    '    '^  '  '''^'       2ur 


where: 


k  =   wave  number. 


c(z)' 
z   =   depth,   and 

P  =   P(r,z,co),   a   function  of   horizontal   range,   depth  and   frequency. 

The  Hankel  transform  allows  one  to  compute  the  changes  of  pressure 
with  a  change  of  depth,  independent  of  the  range.  It  requires  consider- 
ation of   the   vertical    component   of   the  wave  number,    k^: 
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3fp 
3z^ 


+   kzp   =   -TT-   6(z-Zo) 

'^  2lT 


2.4 


where: 
p   =   p(5,rco)    a    function   of    horizontal    wave   number,   horizontal   range   and 
frequency, 


k,^   =  k^  -    e 


^  =  horizontal    wave   number,   and 

k  =   wave  number,    — 7— r. 
c(z) 

This   equation  is   then  solved   by  variation  of   parameters.      Let 


p  =   Wi(z)u(z)   +   W2(z)v(z) 


2.5 


where: 


fz 


Wi(z) 


2tt 


v(z)5(z-Zo)dz 


W. 


uv 


2.5a 


rh 


W2(z)    = 


27T 


u(z)6(z-Zo)dz 


W. 


uv 


2.5b 


'uv 


=   u    -r-  -   v   -7—   (the  Wronskian), 
dz  dz 


2.5c 


and  both  u  and   v  are  solutions   of  the  homogeneous  form   of  equation   2.4. 

Now    if    v(z)   satisfies    the    upper    boundary  condition   and   u(z)   satisfies 
the  lower   boundary  condition,   then: 
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-F  v(Zo) 


w,   =    2-n    W 


uv 


z    >   z, 


z   <    z. 


2.6a 


and 


u(Zo) 


W2   =v 


2^   W 
0 


uv 


I. 


Z     <     Z, 


z   >   z, 


2.6b 


This    means    that    the   transformed   pressure    can   be   written   in   terms    of    depth 
dependent   u(z)   and   v(z)   as: 


P   = 


p      u(Zo)v(z) 
2tt  \l^ 


for   z   <    Zo 


2.7a 


P   = 


p      u(z)v(zo) 
2^        i^W 


for    z   >   z, 


2.7b 


A   convenient   way  of   representing  this  is 


P    =    2^ 


p      u(z>)v(z<) 


W 


uv 


where: 
z>  represents  the  larger  of  z  and  Zq,  and 
z<;  represents  the  smaller  of  z  and  Zq. 

To  obtain  the  actual  pressure  one  must  next  take  the  inverse  Hankel 
transform  of  zero  order: 
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p  = 


pJo(Cr)CdC 


JO 


2.9 


|-00 


2tt 


u(z>)v(z<) 


Wuv 


Jo(Cr)^d5  2.10 


Use  is   made  of   the   fact   that  [Ref.    1] 


Jo(Cr>)   =    1/2 


H^^^(Cr)    +   H^^^(Cr; 


2.11 


and 


H,^2^(^r)e-i^  =   -li^o^\^r) 


2.12 


It  is  then  possible  to  express  the  field  integral  (Equation  2.10)  in 
terms  of  a  Hankel  function  of  the  first  kind  [Ref.  2]: 


P  = 


4^ 


U(Z>)V(Z<)    (1) 


'UV 


2.13 


The  Hankel  function  can  be  approximated  by 


H^"(er) 


V'^F  "' 


<^--^'  (-  I^  *  ^  ) 


2.m 


Restricting  this  , expansion  to  the  first  terra  makes  the  integral  much 
more  manageable  but  restricts  the  minimum  horizontal  range  so  that  E,r 
>>    1  : 


P    = 


-iTT/4 


■7T    V  Trr 


U(Z>)V(Z<)        ifp    _ 

—   e   ^V^dC 


W. 


2.15 


uv 


This  field  integral  represents  the  pressure,  within  the  limits 
imposed  by  truncating  the  expansion  of  the  Hankel  function.  This  integ- 
ral forms  the  basis  for  the  derivation  of  formulae  describing  pressure 
using  two  different  numerical  approaches.  At  this  point,  the  approach  to 
the  normal  mode  solution,  which  will  be  used  in  the  program  El^CT,  takes 
a  different  tack  from  the  approach  which  will  be  used  in  the  FFP  solu- 
tion  and   program. 

B.      NORMAL    MODE    METHOD   FOR   PRESSURE    DETERMINATION 

In  this  paper  the  normal  mode  method  refers  to  the  process  whereby 
approximations  for  the  eigenvalue  of  a  mode  are  successively  refined  by 
numerically  integrating  the  square  of  a  depth  function  across  the  depth. 
A  correction  term  is  developed  from  this  integral  and  is  used  to  refine 
the  approximation   until   the   correction  is   negligible. 

Let   the   guess   at   the  normal    mode   horizontal   wave  number,    E,^,   be   E,. 

For   ^,   v(z)   satisfies 


V"   +   k|  V   =  0 


2.16a 
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where: 


For    5pi,   v^(z)   satisfies 


k|   =   k^  -e 


V   ^   k|nv  =   0 


2.16b 


where 


U2         =      1,2     _p2 


Multiplying  2.1  6a  by  Vj^  and   2.1  6b  by  v  and  subtracting  yields 


vnv"   -   vvn"   +   vvn(k|-k|n)   =   0 


2.17 


Then 


—    (v^v'-vvn')    +   vVn(k^C^-k^+Cn^)    =  0.   O"^ 


dz 


dz 


(v^v'-vv^')    =   (^^C^)   vv^ 


2.18 


Integrating  from   depth   0   to   depth  h: 


(v^v'    -   vv^') 


=  U'  -  C^) 


vv^   dz 


2.19 


or 
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rh 


Vn(h)v'(h)   -    vn'(h)v(h)   -    v^COv'CO)    +   v^COvCO)    =   (C'-^n) 


vv^dz.  2.20 


The  following  argument  holds  for  any  set  of  permissible  boundary 
conditions  but,  in  order  to  clarify  the  discussion,  the  specific  case  of 
pressure  release  at  z  =  0  and  rigid  at  z  =  h  will  be  selected.  The 
function  v(z)  is  not  a  mode,  but  it  can  be  forced  to  satisfy  one  of  the 
boundary  conditions  without  losing  generality.  If  v(z)  ia  defined  to 
satisfy  the  boundary  condition  at  z  =  0,  then,  for  any  permissible  boun- 
dary condition,  the  last  pair  of  terras  on  the  left  hand  side  of  equation 
2.20  is  identically   zero.      Therefore, 


Vn(h)v'(h)   -    Vn'(h)v(h)   =   {^  -C^) 


vvj^   dz'. 


2.21 


For  the   simplest    case   considered  in   this    paper,    the    bottom    boundary   con- 
dition.is   rigid  so  that    V[^'(h)   =   0.      Therefore 


Vn(h)v'(h)   =   ig'' 


a) 


vvj-j   dz 


2.22 


and,   since   v::  v,^   for   E,  close  to   Cn» 


■2  ^r2_i,2   ^    v(h)v'(h: 
/vMz 


AC^  =C^?^   = 
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,2  _   .2       v<h)v'(h) 


2.23 


Besides  proving  useful  in  determining  successively  better  approximations 
of  ^^,  this  same  integral  forms  part  of  the  expression  for  the  contribu- 
tion of  a  mode,  since  it  is  proportional  to  the  derivative  of  the  Wron- 
skian  with  respect   to  wave  number. 

The  Cauchy  integral  enables  one  to  write  equation  2.13  as  the  sum  of 
the  residues  of  all  the  poles  of  the  integrand.  Each  pole  corresponds 
to  a  normal   mode. 


2iT     "— 
n 


•^n^^nHo^   Un'^)Cn 


dW 


2.24 


^n 


where: 

W(Cn)   =   0 
To  derive  an  expression  for   (dW/dC)   multiply  both  sides   by 


Un(h) 


(B  also  equals   ) 


Vnt 


rh 


Un(h)v'(h)  -    Ui!i(h)v(h)   =   ( ^-Cn)(C+5n)  B 


vMz 


2.25 


The   left-hand   side   of   this    equation    gives    the  Wronskian   for    u   and   v, 
and,   as   ^  ^   ^n,   WC^^i)   ^   0,    as   it  should: 
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dW 


2.26 


W(C)-W(Cn) 


fh 


=   B(C+5n) 


v^dz         2.27 


Since    E,    _    Cn»    ^"^^    by    properly   constructing   the    program,    the    constant    6 
can   be  forced  to   be   equal   to    1 , 


dW 

dC 


rh 


=   2^n 


vMz 


2.2S 


and  equation   2.24   can  now    be  written: 


iF     y         ^n^n  (1) 


2.29 


Using  the  exponential  for  the  Hankel  function,  as  in  eqn.  2.13: 


P   =    ^ 

4lT 


y    J_2_     ^nVn       iCnr     "i  tt 


/4 


2.30 


=    _L    y      /     2        ^n^n 

4tt     ^    V^^n"^  -^"^ 


e     "^        i  TT/i4 
e 


'dz 


2.31 
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In    reality,     this    only    represents    the    pressure    in    a    system    for    which 
there  are  no   boundary  losses:    i.e.,   for  self-adjoint   boundary  conditions. 
If    real    losses,    including   scattering,    were   to   be    taken   into   account, 
fh 


the  term 


vMz  must   be   replaced   by  [Ref.    31  = 


where: 


vMz   + 


v^(0)(|^) 


-       ,r2 


v^in){-—:) 


2.32 


A    = 


dv/dz 


U  = 


du/dz 
u 


Ignoring  boundary  losses,  equation  2.31  represents  the  total  pressure 
at  a  point  which  is  at  depth  z  and  horizontal  distance  r  from  a  point  cw 
source  at  depth  Zq.  This  equation  is  the  basis  of  the  computer  program 
EXACT,   which  is   discussed  in   detail   in  Chapter  III. 

C.      DIRECT   EVALUATION    OF   THE    FIELD  INTEGRAL 

The  second  method  of  determining  pressure  at  a  point  requires  calcu- 
lating the  field  integral  itself  and  then  numerically  integrating.  Equa- 
tion 2.15  represents  the  pressure  field  in  terms  of  the  Wronskian  and 
the  source  and  receiver  depth  functions  u(z)  and  v(z).  The  integral  is 
solved  by  replacing  it  with  a  summation  across  an  interval  of  wave 
numbers.  Pressure  at  a  point  is  found  to  be  the  sum  of  the  pressure 
contributions  from   each  horizontal   wave  number: 
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P   =    in^    y-    e  ^L     ^^— —   e        ^      /iAi   AC  2.33 


u(z>),  v(2<),  and  W^y  are  calculated  from  the  known  boundary  conditions 
at  the  surface  and  at  the  bottom.  A^  must  be  small  enough  to  ensure 
accuracy. 

It  is  impractical  to  pursue  this  theme  without  introducing  absorp- 
tion. Evaluation  of  equation  2.3^  when  mA^  =  E,^  produces  W^Jy  =  0  in  the 
denominator,  leading  to  an  infinite  solution.  The  magnitude  of  P^^  then 
depends  on  how  near  mA^  comes  to  E,^.  Introduction  of  absorption  as  an 
imaginary  component  of  horizontal  wave  number  displaces  the  poles  from 
the  real  axis  and  ensures  that  the  Wronskian  does  not  go  to  zero  along 
the  path  of  integration.  This  has  the  effect  of  dispersing  the  energy 
out  of  the  theoretically  limitless  pressure  function  that  results  when 
the  Wronskian  goes  to  zero  at  the  exact  wave  number  corresponding  to  the 
normal    modes. 

The  terms  within  the  summation  sign  of  equation  2.33  have  the  same 
form  as  the  discrete  Fourier  transform,  normally  considered  for  time,  t, 
and  frequency,   oj: 


N-1 

S~     nr    A    ^      i2TTnm/N  ,     _.  ^  _,, 

2_    G(nAco)   e  =   g(mAt).  2.3^ 

n=0 


By  simply  renaming  the   variables,   the  transform   becomes: 
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1_     G(nAC)   e 
n=0 


i2TTnm/N 


=   g(mAr). 


2.35 


This  is  the  basis  for  the  technique  known  as  the  Fast  Field  Program 
[Ref.^].  By  comparing  equation  2.3^  to  equation  2.35,  the  following  can 
be    seen    to   be    true: 


G(nAC)    =    ^    C 


C=nAC 


2.36 


The   pressure   at   range  R    can  then   be   written   as 


Pr   =    PraAr 


iiTT      YtT 


~-    e'^""^^    ACCFFTCG]]. 
mAr 


2.37 


The  capability  of  this  Fast  Field  Program  will  be  limited  by  the  number 
of  range  (and  wave  number)  increments.  The  wave  number  sampling  incre- 
ment, AC,  the  range  resolution,  Ar,  and  the  number  of  points  in  the  FFT 
can   be  seen  to   be  related   by: 


2Trnm 
N 


=    E,r   =    (nAC)(mAr) 


or 


2tt 


=    ACAr, 


2.38 


2.39 
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which  results  from  comparing  equations  2.34  and  2.35.  The  number  of 
sampling  points,  N,  places  practical  limits  on  either  the  maximum  range 
or   the  range   resolution. 

D.      ABSORPTION   AND   BOUNDARY  EFFECTS 

Sound  energy  is  "lost"  by  absorption  in  the  water  and  transmission 
into  the  bottom.  In  both  cases  a  primary  variable  is  frequency  but  the 
angle  of  travel  of  the  mode  and  the  type  of  bottom  must  be  taken  into 
account.  This  paper  limits  discussion  to  a  non-scattering  surface,  water 
volume  and  bottom.  It  also  neglects  energy  that  is  transmitted  into  the 
bottom  and  then  returned  to  the  water  after  refraction  within  the 
bottom. 

Sound  absorption  in  the  water  at  low  frequencies  (under  1  KHz)  is 
due  primarily  to  relaxation  of  the  boric  ion.  Standard  texts  [Refs.  5,6] 
define  absorption  in  this   frequency  range   as: 

a  -    ^^1^  2.  40 

where: 

a  =  absorptive   loss   in  dB/m,   and 

F  =   frequency  in   kilohertz. 
This    loss    can    be    converted    to    the    fractional    loss,    a,    for    a    given    dis- 
tance in  nepers/meter: 

Loss  =   20  log-|o  ^^^   clB 
ar   =   20   logTo  ®"'" 
=   20ar    log-|Q  e 
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=   8.7   or 

a   =   8.7a  dB/meter,   or   a  =    o~^  nepers/meter  2.42 

0.  7 

This    loss   rate    can    be   introduced   into   the   equations    derived    earlier 
by  adding  absorption  as   an   imaginary  component   of   the   wave   number: 

kc   =  k   +   ia  2.il2 

The  exponential   form   provides   the   absorptive   loss   factor: 

ikpR  i(k+ia)R 

e     '-     =   e 


-  e'^"   e""".  2.1.3 


However,    all    of    the    expressions    for    pressure    are    given   in   terras    of    the 

f        r   ■ 
horizontal    range,    r,    not    actual   "path"    length,    R.      Since    -r   =    q,   the   loss 

K  n 

factor  can  be  written  as: 


k 
loss  =  e  "^  ^    ^ .'  2.44 


Bottom  interactions  cause  losses  and  phase  changes  which  must  be 
taken  into  account.  In  their  discussion  of  propagation  loss  using  normal 
modes  ands  an  impedance  boundary  condition,  Koch  and  Lindberg  [Ref.  7] 
begin  their  discussion  by  defining  the  relationship  between  pressure  and 
its   derivative   with  respect   to   depth: 
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!;  =  -z  M  ^ 


where 
P   =    fluid  pressure  in  the   water   at   the   boundary, 

dP 

—   =   derivative   of   pressure  with  respect   to   depth, 

k2   =   vertical   wave  number,   and 

R   =   complex  relfection   (Rayleigh)   coefficient. 

The    Rayleigh    reflection    coefficient    for    the     boundary     between    two 
homogeneous   fluids   is   given  [Ref.    5,    Equation   6.30]  by: 

Pb^b    _    sinBb 

Pw^w         sin9w  ^   ,  , 

R   =    T——  2.46 

Pb^b    ^    sinGb 

Pw^w    ^    sine^ 

where 
c^   =   sound  speed  in  the   bottom   at   the   interface, 
c^  =   sound  speed  in  the  water    at   the  interface, 
Pvj  =   density  of   water, 
Pb   =   density  of   the   bottom, 

9y    =    angle,    measured    from    the    horizontal    at    which    the    wave    strikes 
the 

bottom,   and 
9i-)   =  angle,   measured   from    the   horizontal,   at    which  the   wave   is 
transmitted  into  the   bottom. 
Substituting  this   value  of  R    into  Equation  2.^15  and  simplifying: 
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dP        .,       Pw^w    sine^ 

~r  =   Iky    — — -—   P.  2.i47 


But 


sine^   =   /l-cos^95  2.il8 


Snell's  Law   relates   the  incident   and  transmitted  angles   by: 


°b 
cos    9^,   =    —   cos  8y  2.iJ9 

*^w 


So, 


•iM 


sm   Qb   =   \/^    "    TT   cos^Gw  .•     2.50 

V  w 

The    critical    angle,    e^,    is    defined    as    that    grazing    angle    for    which 
the  transmitted  angle  is    zero.      This   means   that: 


cosew         cos(O)  ,  C5  T 

=    and  —   =    -— ,  2.51 

Cw  Cb  Cw         cosBc' 


where    Qq    is    the    critical    grazing    angle,    measured    from    the    horizontal. 
Therefore,   from   Equation   2.50, 


/■ 


cos^Gu 
^^"'b   -    V/l    -    ^^JT^-  2-52 
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By  defi-nition, 


cos    9,^   =    ^  2.53 


and 

k 


sin    9vv   =    -j^    .  2.54 


Substituting   Equations    2.52  and   2.54  into   Equation   2.4?: 


dP       .,    Pw^w      r       cos^e^ 

dz  PhCh     V  cos^e 


=   ik   ::-^^   \/l    -    ,^^2,      P.  ■  2.55 


In  terms   of   wave  numbers,    Equation   2.53   can   be   used  to   give: 


=   ik    T^^    y1    -    ,.2..^\2.      P  2.56 


dP         .,     PwCw      /  ^2 

dz  PbCb    V  k^GOS^eQ 


or 


dP   =   i  ^^     L2 §! 

dz  PbCb    V  cos^6q 


=   i  -r-^    \/k'  -   — TT7-    P.  2.57 


Equation  2.56  shows  that,  for  any  grazing  angle  less  than  critical, 
the  pressure  derivative  will  be  real  whereas,  for  grazing  angles  greater 
than  critical,  the  derivative  will  be  purely  imaginary.  If  the  pressure 
were   complex,    the    derivative    would   be    complex  and   complicated,   too 
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complicated    to    solve    in   all    cases   with    the   basic   program   design.      Therefore, 
only  the  real   component   of   pressure  was    used  in  the   determination  of   the 
boundary  conditions. 

A  practical  discussion  of  the  physical  significance  of  the  reflection 
coefficient  and  the  relationship  amongst  it,  the  grazing  angle  and  criti- 
cal   angle  is  included  in  Chapter  IV. 


32 


III.       PROGRAMS 

Implementation  of  theory  into  the  computer  programs  required  two 
distinctly  different  methods.  The  program  EXACT  made  accurate  estimates 
of  the  normal  mode  horizontal  wave  numbers  in  order  to  compute  the 
pressure  as  the  sum  of  the  pressure  components  from  each  of  the  modes. 
Absolute  values  of  pressure  are  of  no  concern  in  either  EXACT  or  FFP 
because  the  purpose  is  to  compute  propagation  loss.  It  was  assumed  that 
forcing  function,  F,  had  a  value  of  unity.  The  calculated  pressure  can 
then  be  compared  to  the  pressure  at  one  meter,  and  from  that  ratio,  the 
propagation  loss  is  calculable  using  decibel  ratios.  The  program  FFP  did 
not  compute  horizontal  wave  numbers  of  modes;  it  assumed  a  large  number 
of  wave  numbers  equally  spaced,  calculated  the  pressure  for  each,  and 
summed   the   individual    contributions    to   give   a  final   pressure. 

Both  programs  were  written  using  complex  numbers  and  double  preci- 
sion (64  bit  words),  except  for  subroutine  'eigens'  which  was  restricted 
to  single  precision  because  it  relied  upon  IMSL  routine  EQRTIS  which  was 
available  only  in  single   precision. 

Subroutines  'modesl'  and  'modes2'  provide  many  common  functions  for 
both  programs  EXACT  and  FFP.  They  calculate  the  transmitter  and 
receiver  modal  pressure  functions  by  stepping  upward  or  downward  from 
one  of  the  boundaries,  computing  the  pressure  and  its  derivative  with 
respect  to  depth  by  means  of  a  fourth  order  Runge-Kutta  technique 
[Ref.7].' 
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A.      NORMAL    MODE    PROGRAM 

This  program  is  attached  as  Annex  B.  EXACT  computes  the  approximate 
modal  horizontal  wave  numbers  by  means  of  the  subroutine  'eigens'.  It 
then  refines  the  horizontal  wave  number  estimate  to  the  desired  accuracy 
using  subroutine  'modesi'  or  'modes2'  and  computes  the  pressure  factors 
for  each  mode.  Its  aim  is  to  compute  the  total  pressure  and,  hence,  the 
propagation  loss   at  each  range  of  interest. 

Subroutine  'eigens'  uses  a  matrix  method  [Ref.8]  to  approximate 
all  the  real  modal  horizontal  wave  numbers  as  well  as  the  first  evanes- 
cent mode.  At  the  ranges  of  concern,  it  is  unlikely  that  more  than  the 
first  evanescent  mode  would  be  significant.  The  matrix  of  wave  number 
estimates  is  found  in  a  relatively  short  time  and,  although  not  guar- 
anteed to   be  accurate,  it  is   complete;   no  modes  are  skipped. 

The  next  step  involves  finding  the  exact  values  of  the  horizontal 
wave  number  for  each  mode.  The  procedure  utilizes  the  fact  that,  for 
each  mode,  both  boundary  conditions  for  a  depth  function  (u  or  v  from 
Chapter  II)  will  be  met.  By  setting  the  depth  function  or  its  derivative 
with  respect  to  depth  to  a  known  boundary  condition  at  the  top  or  the 
bottom,  it  is  possible  to  find  a  vertical  wave  number  that  will  result 
in  the  boundary  conditions  being  met  at  the  the  other  boundary.  The 
boundary  with  higher  sound  velocity  will  experience  an  exponential  decay 
in  the  depth  function.  By  starting  at  this  boundary,  it  is  a  relatively 
simple  matter  to  solve  for   vertical    wave  numbers   that  result   in 
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sinusoidal  variations  that  lead  to  solvable  boundary  conditions  as  the 
other  boundary  is  approached.   This  avoids  the  possible  problem  that 
would  exist  if  an  attempt  had  been  made  to  start  at  the  boundary  which 
had  the  lower  sound  velocity.   In  that  case  it  is  possible  that  instead 
of  solving  for  an  exponentially  decaying  function,  the  function  may  be 
exponentially  growing,  and  no  solution  is  possible.   Figure  3.1  illus- 
trates how  a  poor  choice  of  starting  boundary  could  lead  to  an  unsolvable 
situation.   For  each  of  the  modes  in  the  'main'  routine  calls  either 
'modesl'  (when  sound  speed  is  greater  at  the  top  than  at  the  bottom),  or 
'modes2'  otherwise.   It  calls  with  starting  estimates  of  the  modal 
horizontal  wave  number;  the  value  from  'eigens'  the  first  time  and  an 

improved  value  each  subsequent  time. 

'Modesl',    using    the    best    estimate    for    E,^,    computes    the    modal    depth 

function    at    each    increment    of    depth,    starting    at    the    surface    where    the 

depth   function   is   zero   and  its   derivative   with  respect   to  depth  is   set   to 

1.0.      The   integral    (Equation    2.31)    is    calculated   at    each   interval    and   its 

final   value   at   the    bottom   is   used   in  the    correction  term. 

At  the  bottom,  the  depth  function  and  derivative  are  compared  to  the 
boundary  conditions  expected  for  a  mode.  For  a  rigid  bottom,  a  deriva- 
tive of  zero  is  expected.  For  an  impedance  bottom,  the  expected  deriva- 
tive will  be  defined  by  Equation  2.56.  The  error  in  the  derivative,  tog- 
ether with  the  depth  function  and  the  integral  is  used  to  calculate  a 
correction   for   the   modal   horizontal    wave   number. 

If  the  surface  sound  speed  is  greater  at  the  bottom  than  at  the  top, 
subroutine  'modes2'  is  called.  It  performs  the  same  function  as  'modes!' 
but    starts    at    the    bottom,    where    the    depth    function    and    derivative    are 
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Figure  3.1.   Possible  SoluCions  to  Depth  Function  Starting  at: 
a:  Boundary  with  Exponential  (Convergent)  and 
b:  Boundary  with  Sinusoidal  Changes  (Possibly  Divergent) 
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set  to  meet  the  boundary  conditions  for  an  impedance  bottom  (Equation 
2.56).  A  correction  is  made  after  the  subroutine  has  stepped  to  the 
surface  and  comparison  has  been  made  with  modal  boundary  condition 
expectations;    in  this    case   depth   pressure   function   equal   to   zero. 

'Modesl'  or  'modes2'  will  be  repeatedly  called  until  E,^  has  been 
satisfactorily  estimated.  The  'main'  routine  then  computes  the  'incom- 
plete' modal  pressure  from  the  modal  depth  functions,  the  integral,  and 
the  horizontal  wave  number,  E,^,  (Equation  2.31)  for  each  mode,  storing 
them  for  later  use.  It  is  termed  'incomplete'  because  there  is  no  range 
dependency  at   this  point. 

After  the  last  mode  has  been  processed  by  'modes!'  or  'modes2', 
range  is  introduced  and  partial  pressures  for  each  mode  are  calculated, 
keeping  an  updated  sum  (complex)  of  all  the  modal  pressure  factors. 
These  pressure  factors  are  then  converted  to  propagation  loss.  It  is  a 
simple  matter,  once  all  the  eigenvalues  have  been  found,  to  compute  the 
propagation  loss  for  a  series  of  ranges,  as  can  be  seen  in  Appendix 
B,    Program  EXACT. 

B.      FAST   FIELD   PROGRAM 

This  program  was  written  with  the  intent  of  completing  a  Fast 
Fourier  Program,  FFP,  which  would  utilize  a  Fast  Fourier  Transform,  FFT, 
subroutine.  However,  the  program  was  only  completed  to  the  point  where 
pressure  factors  were  computed  for  individual  ranges  (as  in  EXACT), 
rather  for  a  large  number  of  ranges  (as  is  the  intent  of  an  FFP).  The 
computer    program   FFP,    although  incomplete,   is  included  as  Appendix  C. 
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It  was  mentioned  previously  that  EXACT  and  FFP  are  provided  many 
common  functions  by  'modesl'  and  'modes2'.  There  are  also  some  basic 
differences  in  the  way  the  subroutines  are  used  by  the  two  programs.  In 
FFT,  since  no  exact  wave  numbers  are  to  be  calculated,  there  is  no 
requirement  for  subroutine  'eigens'.  For  the  same  reason,  there  is  no 
requirement  to  find  corrections  to  the  wave  numbers,  so  'modesl'  and 
'raodes2'  are  used  differently.  In  fact,  both  subroutines  are  required  by 
FFP. 

Using  the  incremental  wave  number,  nA^,  subroutine  'modesl'  starts 
at  the  surface  with  the  same  boundary  conditions  as  in  EXACT  and  steps 
down,  computing  the  pressure  function  and  the  derivative,  through  the 
depth  of  the  upper  of  the  transmitter  or  the  receiver,  and  stops  at  the 
lower  of  the  two.  Values  are  required  at  the  lower  level  for  later  use 
with  the  values  from  'modes2'  to  calculate  the  Wronskian.  Then  'modes2' 
starts  at  the  bottom  with  the  same  boundary  conditions  as  EXACT  and 
steps  upward  until  the  pressure  function  and  derivative  have  been  calcu- 
lated  for   the   lower   level. 

The  'main'  routine  then  uses  the  pressure  functions  and  derivatives 
to  calculate  a  partial  pressure  (Equation  2.33)  for  that  wave  number  in- 
crement. Each  partial  pressure  is  summed  (integrated)  so  that,  when 
the  last  wave  number  has  been  processed,  the  result  is  the  total  pres- 
sure. 

As  in  EXACT,  FFP  computes  the  propagation  loss  for  each  range. 
Implementation  of  the  FFT  would  make  it  unnecessary  to  step  through 
ranges  as  must  be  done  with  EXACT;  propagation  loss  would  be  calculated 
for   as  many  ranges   as   there   are   FFT   points. 
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IV.    TESTS    RUN    ON    THE    PROGRAM 

Because  it  was  difficult  to  gauge  the  correctness  of  the  program  or, 
more  to  the  point,  the  veracity  of  this  author's  ideas,  it  was  considered 
prudent  to  check  the  results  at  each  stage  of  development.  Although  the 
process  was  awkward  and  time-consuming  it  proved  necessary  and  resulted 
in  some  unexpected  benefits;  the  results  of  some  of  the  tests  provided 
graphical  insight  into  the  phenomenon  of  sound  propagation.  By  starting 
with  isospeed  conditions,  verifications  were  simplified.  Solutions  for 
horizontal  wave  numbers  which  resulted  from  the  computer  programs  were 
checked  against  'correct'  solutions  from  simple  formulae  for  the  non- 
absorption,  rigid  bottom  system,  as  well  as  for  the  system  that  included 
absorption  in  the  water.  Once  the  impedance  bottom  was  introduced,  an 
analytic  approach  was  required  to  decide  if  corrections  were  applied 
properly.  A  simple  expedient  to  check  on  the  program  at  any  point  in- 
volved plotting  propagation  loss  against  horizontal  range  for  pairs  of 
modes  and  observing  the  interference  distance  between  nulls.  This  dis- 
tance, when  compared  to  the  theoretical  distance  would  highlight  an 
error  if  there  were  an  anomaly.  Another  check  involved  the  observation 
that,  for  a  shallow  channel  such  as  that  used  in  the  model,  losses  are 
generally  spherical  at  short  ranges  and  cylindrical  for  far  ranges.  The 
smoothed  propagation  loss  curves  could  be  expected  to  lie  roughly  along 
curves   predicted  on  this   basis. 

The  first  tests  involved  studying  the  basic  building  blocks  of  the 
program.      As   the   data  were  read  from  input   files,    they   were    printed   into 
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a  new  file  to  ensure  that  correct  information  was  being  utilized. 
Another  simple  check  involved  printing  a  matrix  of  the  wave  numbers  as 
they  were  computed.  To  illustrate  the  importance  of  such  a  basic  check, 
it  is  worth  remarking  that  the  half-increment  method  required  that  the 
wave  number  be  calculated  for  each  half-increment  of  depth,  including  a 
depth  of  zero.  Since  the  variation  of  Fortran  used  does  not  support  a 
matrix  with  a  zeroth  element,  a  printout  proved  very  useful  in  visualiz- 
ing the  situation  and  pin-pointing  a  subtle  error.  In  many  cases  it  was 
only  by  printing  out  all  of  the  variables  after  each  step  that  errors 
could   be  identified  and  remedied. 

A.      RIGID   BOTTOM   WITH  NO  ABSORPTION 

Subroutine  'eigens'  was  first  checked  for  correctness  by  comparing 
its  constant  speed  solution,  for  each  mode,  to  a  solution  which  was 
known  to  be  accurate.  This  accurate  solution  was  based  on  the  concept 
that,  for  a  constant  sound  speed,  with  rigid  bottom  and  no  absorption, 
the  pressure  function  would  be  sinusoidal,  satisfying  both  boundary  con- 
ditions. This  means  that  the  horizontal  wave  number  for  each  mode  (n 
always  odd)  can  be  defined  in  terms  of  the  wave  number,  k,  and  the  water 
depth,   H: 


cA  =  k^  -  (^y-  ^-T 


Because  the  subroutine  could  only  be  run  in  single  precision,  high  reso- 
lution was  not  expected.  Table  4.1  shows  that  resolution  increased  lin- 
early  with  the  number   of   depth  increments.      For   example,    for   mode   2,   the 
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TABLE  4.1.   SUBROUTINE  'EIGENS'  -  COMPUTED  SQUARES  OF  WAVE  NUMBERS  AND 

ERRORS  RESULTING  FROM  VARIATION  OF  NU^ffiER  OF  DEPTH  INCREMENTS 


MCDE  NUMBER 
2 


NUM  OF 
INCRMT 


VALUE 
ERROR 


VALUE 
ERROR 


VALUE 
ERROR 


VALUE 
ERROR 


25     0.2144624294 
-0.0000906215 


0.1959892982    0.1525131932 
-0.0009006298   -0 . 00328778S3 


0.0237859429 
-0.05145329C6 


50     0.2144174816 

-0.C000456736 


0. 1955415049 
-0.0C0452S355 


0. 1508744489 
-0.0016490440 


-0.0209414381 
-0.0157259096 


100    0.2143947361 
-0.0000229281 


0. 1953157087 
-0.0002270403 


0. 1500510807 
-0.0008256757 


-0.0305258121 
-0.0071405356 


200    0.2143832949 
-0.0000114870 


0.1952023429    0.1496385139 
-0.0001136746   -0.0004131089 


-0.0342955159 

-0.003371S318 


400    0.2143775572 
-0.0000057492 


0.  1951455441 
-0.0000568758 


0. 1494320249 
-0.0002066200 


•0.0360238122 
-0.0015435355 


800    0.2143746840 
-0.0000028760 


0.  1951171158 
-0.0000284475 


0. 1493287310 
-0 . 0001033261 


-0.0368554579 
-0.000811839S 


1600   0.2143732463 

-0.0000014384 


0.1951028945 
-0.0000142262 


0. 1492770720 

-0.0000516671 


-0.0372637914 
-0.0004035563 


3200   0.2143725272 
-0.0000007193 


0.  1950957821 
-0.0000071137 


0. 1492512396 
-0.0000258346 


•0.0374661565 
-0.0002011912 


E>'J\CT    0.2143718079         0.1950386633         0.1492254049       -0.03  76673477 
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error  using  25  increments  was  .000900,  and  the  error  was  halved  by  in- 
creasing the  number  of  increments  to  50.  This  trend  continued  through 
the  run  that  used  3200  increments  and  resulted  in  an  error  of  0.000007. 
It  was  decided  that  an  error  of  0.00002,  which  resulted  from  the  use  of 
100  increments,  was  acceptable.  The  value  of  each  horizontal  wave 
number  determined  in  subroutine  'eigens'  was  intended  only  as  a  starting 
point  for  the  more  accurate  subroutines  'modes'!'  and  'modes2'.  The 
accuracy,  therefore,  has  significance  only  in  that,  if  two  modes  are 
closer  together  than  single  precision  accuracy  limitations,  there  would 
be   a   chance  that   a  mode   will   be   missed   completely. 

Subroutine  'eigens'  was  also  checked  using  sound  speed  profiles  that 
varied  linearly,  both  positively  and  negatively,  with  depth.  Since  no 
clear-cut  means  of  assessing  the  error  was  available,  no  conclusions 
could  be  drawn,  and  results  are  published  for  only  the  positive  gradient 
(Table  i4.2a).  In  addition,  by  using  a  sound  speed  profile  that  produced 
two  strong  ducts,  it  was  possible  to  show  that  'eigens'  was  capable  of 
distinguishing  between  two  horizontal  wave  numbers  which  are  very  close 
together  (Table  4.2b).  The  solutions  are  not  precise  but  no  modes  are 
missed  and  the  necessary  information  is  provided  to  the  subsequent  sub- 
routines  where  the   values   are   calculated  to  the  required   precision. 

Next,  normalized  modal  pressures  at  a  large  number  of  depth  incre- 
ments were  calculated  using  'modes2'.  These  data  were  combined  in 
Figure  4.1  and  show,  for  each  mode,  the  relative  amplitude  of  the  modal 
pressure  function.  The  profiles  for  all  three  real  modes  as  well  as  the 
first  evanescent  mode  are  plotted  together  and  care  must  be  taken  to 
note  that   they  represent   the  modal   characteristics,    not    the   actual 
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TABLE    4.2.       SUBROUTINE    'EIGENS'    -    SQUARES    OF    MODAL    HORIZONTAL 

WAVE    NUMBERS    FOR      A.    POSITIVE   GRADIENT,    AND    B.    DUAL 
DUCT    FOR    50   HZ   SOURCE    AND    500   METERS    DEPTH. 


MODE  POSITIVE  .DOUBLE 

NUM  GRADIENT  DUCT 


1 

0 

.2023466395 

2 

0 

.2017806065 

3 

0 

.2015439681 

4 

0 

.2009845610 

5 

0 

.2001803173 

6 

0 

.1991866971 

7 

0 

.1980052050 

8 

0 

.1966276551 

9 

0, 

.1950384284 

10 

0, 

.1932226702 

11 

0, 

.1911713327 

12 

0. 

.1888803823 

13 

0, 

.  1863476142 

14 

0. 

, 1835693580 

15 

0. 

,1805375878 

16 

0.. 

.  17723S1208 

17 

0. 

,1736509107 

18 

0, 

,  1697523042 

19 

0, 

,1655174308 

20 

0. 

,1609207800 

21 

0. 

1559344080 

22 

0. 

,1505243249 

23 

0. 

,1446458558 

24 

0. 

1382386306 

25 

0. 

1312213613 

26 

0. 

1234851760 

27 

0. 

1148818518 

28 

0. 

1051995705 

29 

0. 

0941103985 

30 

0. 

0810451735 

31 

0. 

0648246885 

32 

0. 

0419336933 

33 

-0. 

0276373313 

34 

0 

0.2094159534 
0.20922753S2 
0.20S35015S0 
0.20S2S2S212 
0.2075239702 
0.20£5714:'50 
0.2054225531 
0.2C40741S97 
0 . 2C252209S5 
0.2007616530 
0. 19S7S73347 
0. 1965925365 
0. 1941702449 
0.  1915113532 
0. 1S36060596 
0. 1S5442"397 
0. 1320079592 
0. 1732360324 
0. 1742535647 
0. 1699C3S333 
0.  1551959531 
0. 1501C33C95 
0. 1545394076 
0. 1435057637 
0. 142093  6161 
0. 1349764974 
0.  1271523591 
0. 1184323554 
0. 1037542301 
0. 0976346549 
0.0347116344 
0. 0637321079 
0.0469792092 
•0. 0193930514 
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pressure   contributions.      The   figure   shows    that,    for   rigid   bottom,    no   absorp- 
tion,   iso-speed    conditions,    the    profiles    meet    the    boundary    conditions    at 
the    top,    where   pressure   is    zero,    and  at   the   bottom,   where  the   derivative 
is   zero.      The  same   test    procedure   was    repeated    using  subroutine    'modesl' 
and  identical   results    were  obtained. 

To  illustrate  the  relative  degree  to  which  each  mode  is  stimulated, 
depending  on  the  depth  of  receiver  and  transmitter,  a  trial  run  was  made 
with  a  transmitter  at  depth  13  meters  and  receiver  at  various  depths. 
Figure  4.2  depicts  the  relative  strength  of  each  mode  at  each  depth. 
The  second  mode  is  most  strongly  stimulated  and  mode  1  is  least 
strongly  stimulated,    as    could   be   anticipated   from  Figure    4.1. 

Routine  FFP  does  not  use  exact  normal  mode  horizontal  wave  numbers, 
but  boundary  conditions  must  remain  consistent  and  sinusoidal  variations 
still  occur  at  a  rate  defined  by  the  vertical  wave  number,  k^^,  which,  in 
turn,  is  a  function  of  the  horizontal  wave  number,  ^.  The  results  of 
subroutine  'modesl'  were  compared  to  the  values  of  sin(l<2Z>)  and  the 
results  of  'modes2'  were  compared  to  the  values  of  cos(k2(H-z<),  where  H 
is  the  water  depth  and  z>  and  z<  are  the  depths  of  the  receiver  and 
transmitter.  The  comparisons  showed  that  the  subroutines  worked 
properly. 

At  this  point,  the  solution  using  calculated  horizontal  wave  numbers 
in  EXACT  was  examined  for  accuracy.  It  was  decided  that  the  accuracy 
criterion  would  be  based  upon  the  worst  case  which  allowed  for  all 
errors  to  be  cumulative.  It  can  be  seen  in  Equation  2.31  that  an  error 
in  ^^  will  exhibit  itself  directly  in  the  square  root  and  in  the  exponen- 
tial.     A   small   error    causes   only  a  small   error  in  the   square  root    but   it 
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Figure  4.1.   Modal  Pressure  Function  for  Three  Real  Modes 
and  Evanscent  Mode,  Individually  Normalized. 
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RELATIVE  PRESSURE  CONTRIBUTION 


Figure  4.2.   Pressure  Contributions  Normalized  to  the  Maximum  Value 

of  the  Most  Strongly  Excited  Mode  for  Source  at  13  Meters 
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is  critical  in  the  phase  component.  It  was  categorically  decided  that 
the  maximum  allowable  phase  error  for  each  mode  would  be  0.1  radians, 
and  that  the  maximum  range  would  be  100,000  meters.  This  limits  the 
maximum   error   in   ^n   ^^    "lO"^   ^O'"   ^H   wave   numbers.. 

In  assessing  possible  real  errors,  a  vector  plot  was  made  for  each 
of  several  solutions.  A  constant  adjustment  was  then  introduced  to  each 
of  the  horizontal  wave  numbers,  giving  them  an  artificial  error.  It  can 
be  seen  in  Figure  4.3  that  for  short  ranges,  a  large  error  is  admissible 
and,  even  at  -50,000  meters,  the  limit  imposed  by  the  accuracy  criterion 
is   well    within  the   allowable   error. 

Vector  plots  were  also  used  to  illustrate  how  a  small  change  of 
range  can  cause  a  dramatic  change  in  pressure.  Figure  ^.Ua  shows  the 
change  of  propagation  loss  from  830  to  840  meters.  Most  of  the  change 
is  involved  with  phase  and  the  amplitude  change  is  small.  This  can  be 
seen  best  in  Figure  4.4b  where  the  range  change  is  only  one  meter.  It 
should  be  noted  too,  that  these  errors  have  been  made  cumulative 
whereas  it  is  unlikely  that  the  errors  would  ever  be  that  way.  One 
would  therefore  expect  the  total  error  to  be  less  than  that  illus- 
trated. 

A  plot  of  loss  with  range  for  a  rigid  bottom,  Figure  4.5,  showed  a 
general  20  log  R  increase  in  propagation  loss  close  to  the  source  and  a 
change  of  10  log  R  as  the  distance  increased.  Superimposed  on  this  is  a 
strong  interference  pattern  from  the  three  real  modes.  Sudden  changes 
of  pressure  with  changing  range,  as  were  noted  in  Figure  4.4,  are  seen  as 
rapid  fluctuations  or  propagation  loss  in  this  plot.  A  decisive  check  on 
the    general    reliablity    of    the    program    to    this    point    involved    excluding 
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Figure  4.3.   Pressure  Components  for  Three  Real  Modes  I llus tracing 
Cumulative  Error  from  Errors  in  Horizontal  Wave  Number 
at  50  Meters  (Upper)  and  50  Km  (Lower). 
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Figure  A. 5.   Propagation  Loss  for  Source  at  10  Meter  Depth 
and  Receiver  at  37  Meters  in  50  Meter  Ocean. 
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all  but  two  real  modes  and  plotting  the  resulting  propagation  loss 
curves.  For  each  pair  of  modes,  the  interference  distance  can  be  shown 
to   be  related  to  the   horizontal    wave   numbers   by: 

^m    <in 
where: 

R    is   the   horizontal   distance   between  nulls, 

and   E,^   and   E,^  are  the   real   horizontal    wave   numbers. 

Comparison    was    made     between    the    calculated    interference    distance    and 

the   observed   distance    for    all    three   interference   combinations   and   proved 

to    be    correct   in   all    cases.      Figure    4.6   shows    the   result    when    the    first 

and   third   modes    were   used.      The   cycle   distance   was    calculated   to   be    96.5 

meters,    which  agrees    with  the   distance  from   the   plot. 

B.      RIGID   BOTTOM   WITH   ABSORPTION 

Absorption  due  to  water  (low  frequency),  although  normally  very 
small,  was  made  artificially  high  to  check  on  the  mechanics  of  the 
program.  With      absorption      set      to      an      arbitrary      value      of      0.0005 

nepers/meter,  it  was  noted  that  the  imaginary  component  of  the  wave 
number  increased  with  the  increasing  grazing  angle  associated  with  higher 
modes:  0.000505  for  mode  1,  0.000555  for  mode  2,  and  0.000727  for  mode 
3. 

To  ensure  that  output  losses  were  consistent  with  those  input,  other 
propagation  loss  curves  were  drawn  using  a  variety  of  rates  of  absorp- 
tion. Figure  4.7  shows  that  the  difference  in  loss  at  any  given  range  is 
8.7  times  the  difference  of  absorption  rates  (as  expected  from  Equation 
2.42).        For     example,     a     change     from     0.0000001     nepers/meter    to    0.0001 
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Figure  4.6.   Interference  Pattern  Between  Modes  1  and  3 
for  Source  at  10  Meters  and  Receiver  at 
37  Meters  in  50  Meters  of  Water  for  50  Hz. 
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Figure  4.7.   Propagation  Loss  for  Absorption  Kate  0.000001  Nepers/Meter 
(Upper)  and  0.0001  Nepers/Meter  (Lower). 
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Figure  4.8.   Propagation  Loss  as  a  Function  of  Depth  for  Sources  a; 
a.  1  Meter,  b.  5  Meters,  c.  15  Meters,  d.  25  Meters, 
e.  40  Meters,  and  f.  45  Meters. 
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nepers/meter  causes  a  change  of  propagation  loss,  at  10  kilometers,  of 
9.0  dB.      One  would   have   expected  a   change  of    1.0  nepers   or    8.7  dB. 

A  series  of  runs  was  made  to  illustrate  how  pressure  varies  with 
receiver  depth  at  a  set  range  for  a  variety  of  source  depths.  From  the 
graphs  (Figure  4.8)  it  is  important  to  note  that  for  a  source  or  a 
receiver  at  the  surface,  the  pressure  will  be  zero '  and  the  propagation 
loss  infinite  because  none  of  the  modes  is  stimulated.  This  does  not 
change  when  an  impedance  bottom  is  introduced,  but  it  would  change  if  a 
real   surface   were   considered. 

Using  these  same  depth  profiles,  it  was  possible  to  further  verify 
the  program  by  checking  for  reciprocity.  The  propagation  loss  is,  for 
example,  the  same  for  a  combination  of  a  source  at  5  meters  and  receiver 
at   25   meters   as   for   a  source   at   25  meters   and  receiver    at    5   meters. 

The  final  set  of  tests  for  a  rigid  bottom  involved  the  use  of  the 
FFP.  At  the  point  of  testing,  the  Fourier  transform  (FFT)  had  not  yet 
been  introduced  so  pressure  factors  were  being  calculated  for  individual 
ranges.  Figure  4.9  shows  that,  for  a  zero  absorption  loss,  the  only  hor- 
izontal wave  numbers  that  contribute  to  the  total  are  at  the  exact  mode 
values.  Small  changes  of  the  horizontal  wave  number  cause  serious 
changes   in  the    value   close  to  the   exact   mode   value. 

Introduction  of  absorption  causes  the  energy  to  be  dispersed  across 
the  spectrum.  There  are  small  contributions  from  across  the  wave 
number  spectrum  and  flattening  at  the  modes.  Higher  absorption  causes 
greater  dispersion  of  the  'pressure'.  The  integral,  with  absorption, 
can  be  seen  to  fluctuate  in  Figure  4.10,  staying  close  to  zero  until  the 
first    mode     is     approached.  The     integral     increases    quickly,    but     not 
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abruptly  when  a  mode  is  reached.  A  sharp  decline  is  notable  when  the 
1 80  degree  phase  change  occurs  at  the  mode.  The  final  value  represents 
the  final   pressure  at   the  receiver. 

Various  checks  were  made  to  ensure  proper  functioning  of  the  FFP 
program  for  a  given  range.  By  changing  the  number  of  wave  number  incre- 
ments, it  was  found  that  the  maximum  amplitude  for  each  mode  changed 
(Figure  4.11)  but  that  the  integrated  pressure  was  identical  when  the 
number  of  increments  was  changed  from  1024  to  2048.  Pressures  agreed 
with  one  another  to  the  seventh  decimal  and  fifth  significant  digit. 
Given  more  time,  it  would  have  been  interesting  to  see  how  few  samples 
could  be  taken,  for  a  single  range,  before  significant  errors  were  intro- 
duced. 

Figure  4.12  shows  the  effect  of  changing  absorption  rates.  However, 
the  change  of  pressure  at  a  given  range  does  not  correlate  with  the 
change  of  absorption  rate.  This  indicated  that  the  FFP  routine  was  in- 
correct in  some  way  and  further    checks   were  required. 

The  integrated  pressures  from  this  program  were  then  compared  to 
the  values  obtained  from  the  program  EXACT  for  different  absorption 
rates  and  ranges.  As  can  be  seen  in  Table  4,3,  agreement  was  poor. 
Because  of  time  constraints,  it  was  decided  to  abandon  the  FFP  and  con- 
centrate on  solutions    using  the  EXACT   method. 

MOTE:  A  great  deal  of  time  had  been  spent  on  FFP  trying  to  find  a  solu- 
tion for  a  fixed  range,  a  prerequisite  to  implementation  of  the  'FFP' 
subroutine.  At  the  time  that  work  with  the  FFT  was  first  suspended,  it 
gave  inconsistent  results  that  did  not  meet  any  of  the  testing  criteria. 
A    small    change    of    increment    size,    from    2047    to    2048   increments,    caused 


TABLE    4. 3.      CALCULATED   PRESSURE    FOR    PROGRAM   FFP    AND   EXACT. 


RANGE 

ABSORPTION 

(METERS) 

(NEPERS/METER) 

1000 

.   0.0000001 

1000 

0.00001 

2000 

0.0000001 

2000 

0.00001 

PRESSURE  FACTOR 

FFP  EXACT 

0.00172867  0.00139872 

0.00128178  0.00138953 

0.00001486  0.00139504 

0.00010503  0.00138101 
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outlandish  changes  in  the  integrated  pressure.  The  changes  of  pressure 
resulting  from  a  changed  absorption  rate  did  not  correspond,  even  rem- 
otely, with  the  changes  expected.  At  that  point,  priorities  were  shifted 
and  it  was  decided  to  concentrate  on  solving  the  'EXACT'  program  prob- 
lems at  the  expense  of  the  'FFP'.  At  a  later  date,  too  late  to  be  prac- 
tically helpful,  some  of  the  FFP  problems  were  rectified  and  the  above 
mentioned  inconsistencies  were  remedied.  It'  was  possible  to  relate 
propagation  loss  to  absorption  rate  and  it  was  shown  that  the  number  of 
samples  has  a  negligible  effect  on  the  final  pressure.  However,  the 
discrepancy  between  the  EXACT  solution  and  the  FFP  solutions  still  exist 
and  the  FFP    is  incomplete. 
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C.      IMPEDANCE    BOTTOM 

The  introduction  of.  an  impedance  bottom  to  replace  the  perfectly 
reflecting  rigid  bottom  made  it  possible  to  better  predict  the  propaga- 
tion loss  experienced  in  reality.  This  was  done  by  defining  the  bottom 
in  terms  of  its  density  and  the  sound  speed  in  the  bottom  at  the  inter- 
face. After  continued  failure  of  the  program  to  solve  for  modes  whose 
grazing  angles  were  far  below  the  critical  angle  and  the  eventual  reso- 
lution of  that  problem,  a  series  of  tests  was  devised  to  try  to  optimize 
the  program.  These  tests  also  provided  a  convenient  means  of  illustrat- 
ing how  various  angles  of  incidence  associated  with  the  different  modes 
resulted  in  the  varying  phase  changes  and  amplitudes  of  reflected 
energy. 

As  with  the  rigid  bottom  case,  it  was  possible  to  plot  the  modal 
pressure  as  a  function  of  depth,  from  the  reflecting  surface  to  the  imp- 
edance bottom.  Figure  4.13  shows  how  the  pressure  reaches  a  maximum 
well  above  the  bottom  for  mode  1.  The  grazing  angle,  in  this  case,  is 
less  than  the  critical  angle  and  illustrates  the  phase  advance  at  the 
bottom.  The  test  was  repeated  for  strong  positive  and  negative  gradi- 
ents  but   showed  little   change. 

The  first  set  of  tests  to  check  on  final  accuracy  involved  using  a 
critical  angle  of  0.3  radians  and  was  intended  to  find  the  optimal 
number  of  increments.  It  had  been  determined  earlier  that  accuracy  to 
the  ninth  decimal,  or  better,  was  obtainable  (iso-speed)  for  all  modes 
using  200  depth  incements/mode.  At  100,  the  accuracy  dropped  to  the 
seventh    decimal    and    at    50,    answers    were    accurate    to    only    the    fourth 
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Figure  4.13 


RelaCive  Pressure  Function  for  3  Real  Modes  for  Critical 
Angle  0.3  Radians,  Isospeed  Water  Conditions,  50  Hz  Source 
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decimal.  Table  4.4  shows  horizontal  wave  numbers  that  were  calculated 
using  100,  200,  400  and  600  depth  increments/mode  in  combination  with 
correction  factors  of  1,  2,  5  and  10.  The  correction  factor  was  divided 
into  the  mode  wave  number  correction  term  of  subroutines  'modesi'  and 
'modes2'to  limit  the  size  correction  and  thus  prevent  instability  . 
These  correction  factors  were  found  necessary  for  modes  whose  grazing 
angles  were  less  than  the  critical  angle  (in  this  case  only  mode  1). 
Without  a  correction  factor,  the  program  would  fail  to  converge  to  an 
acceptable  solution   but    would   become   divergent. 

A  major  problem  in  deciding  the  best  choice  of  variables  came  about 
because  there  was  no  correct  answer  against  which  to  appraise  the 
values.  It  was  reasoned  that  800  increments  per  mode  and  a  correction 
factor  of  10  would  give  the  most  accurate  answer  and  the  problem  nar- 
rowed down  to  finding  a  combination  that  gave  an  acceptable  accuracy 
with  a  reasonable  number  of  iterations.  Having  decided  that  200  incre- 
ments/mode was  sufficient  and  necessary  for  the  rigid  bottom  case,  it 
became  a  question  of  how  large  the  correction  factor  could  be.  The 
choice  was  further  narrowed  to  a  correction  factor  of  2  with  an  error  of 
9.8  X  10~®  and  40  iterations  required  and  a  factor  of  5  with  an  error  of 
3  X  10~^  and  89  iterations.  Because  the  error  was  borderline  for  the 
correction  factor  2,  the  final  choice  was  200  intervals/mode  with  a  cor- 
rection factor   of    5. 

Using  the  inputs  as  decided  above,  an  effort  was  made  to  correlate 
actual  accuracy  with  that  specified  as  a  requirement.  Results  for 
trials    on   all    modes    using  accuracy   criteria   from    10~^   to    10"^^  follow    in 
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TABLE  4.4 


PROGRAM  'EX^XCT'  -  REAL  AND  IMAGINARY  COMPONENTS  OF  -C^ 
FOR  MODES  1  TO  4  RESULTING  FROM  VARIATION  OF  NUMBER 
OF  DEPTH  INCREMENTS  AND  'CORRECTION'  FACTORS. 


COR  i:;c  ITER 

1   100   90 


0.212512534 
0.000450164 


0. 195507949 
-0.003339831 


0. 150C35435 
•0.010172100 


C.002S57414 
0.037779947 


100   38 


0.212512955 
0.000450172 


0. 195507949 
-0.003339331 


0.  150035436 
-0.010172100 


0.002567414 
0.037779947 


100   39 


0.212513034 
0.000450172 


0. 195507949 
-0.003339331 


C  .  150035433 
-0.-010172100 


0.C02357414 
0.037779947 


10  100  135 


0.212513060 
0.000450170 


0. 195507949 
-0.003339331 


0.  15C035436 

-0.C1Q1721C0 


0.002367414 
0.0377-9947 


200   74 


0.2125123S5 
0.000450155 


0. 195507947 
-0.003339330 


0. 150035443 

-0.010172094 


0.002S67435 
C. 037779560 


200   40 


0.212512963 
0.000450171 


0. 195507947 
-0.C03339330 


0.  1SC035443 
■C.010172094 


;02S57435 
;377-9560 


200   S9 


0.21251303: 
0.00045017; 


0. 195507947 
-0.003239830 


0. 15003544: 
-0.01017209.: 


3.00255-435 
0.03'7779560 


10  200  195 


1   400   67 


400 


400   89 


10  400  1=6 


800   72 


800   39 


800   90 


10  800  187 


E.XACT 


0.212513063 
0.000450170 

0.212512335 
0.0C0450155 

0.212512969 
0.000450171 

0.212513033 
0.0004501-2 

0.212513055 
0.000450170 

0.212512835 
0.000450165 

0.212512959 

0.000450171 

0.212513039 
0.000450171 

0.212513066 
0.000450170 

0.214371320 


(RIGID  BOTT)   0.00050  5  3  40 


0. 195507947 
-0.003339330 

0. 195507945 
-0.003339330 

0. 195507946 

-0.003339830 

0. 195507946 
-0.003339330 

0.195507945 
-0.003339330 

0. 195507946 
-0.003339830 

0.  195507945 
-0.003339830 

0.195507946 
-0.003339830 

0.195507946 
-0.003339830 

0.195038818 
0.000555289 


0.  150035443 
-0.010172094 

0.  150035443 
-0.010172093 

0.  150035443 

-0.C10172C93 

0.  15C035443 
-0.010172093 

0. 150035443 
-0.010172093 

0. 150035443 
-0.010172C92 

0.  150035443 
-0.010172092 

0. 150035443 
-0. 010172092 

0.  150035443 
-0.010172092 

0. 149226333 
0.000725949 


0.C02S67435 
0.037779650 

0.002857437 
C. 037779642 

0.CO2357437 
0 . 03  7779642 

0.002367437 

C. 002367437 
0.037779542 

0.002S67437 
0. 037779641 

0. 002857437 
0.037779541 

0.002867437 
0.037779641 

0.002367437 
0.037779541 

0.002867437 
0.037779640 
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Table  4.5.  Assuming  that  the  answers  obtained  using  10"^^  had  the  least 
error  and  could  then  serve  as  a  standard,  errors  were  calculated  as  the 
difference  between  the  standard  and  each  answer.  It  was  found,  in  all 
cases,  that  the  actual  errors  were  slightly  less  than  those  specified  by 
the  accuracy  criteria.  This  was  to  be  expected  since  corrections  were 
made  to  make  the   error  less   than  the   designated  amount. 

The  final  check  involved  the  problems  associated  with  low  modes 
whose  grazing  angles  were  less  than  the  critical  angle.  A  series  of 
calculations  was  made  to  find  the  number  of  iterations  required  for 
different  critical  angles  (different  sound  speed  ratios)  using  a  100 
meter  deep,  isovelocity  ocean  that  supported  seven  real  modes.  It  was 
also  intended  to  uncover  further  problems  associated  with  large  critical 
angles.  At  the  same  time,  the  squared  values  of  horizontal  wave  numbers 
for  the  seven  modes  were  calculated  for  each  critical  angle  and  are 
shown  in  Table  4.6.  For  angles  close  to  the  grazing  angle  for  the  first 
mode,  0.1  radians,  the  program  failed  using  both  subroutines  'modesi'  and 
'modes2'  during  these  runs.  This  was  unusual  because  previous  data  had 
been   collected  for   almost  identical    conditions. 

The  fraction  of  energy  reflected  and  the  phase  change  represented  by 
the  Rayleigh  coefficient,  were  calculated  for  each  mode  for  each  criti- 
cal angle  and  are  shown  in  Figure  4,14  and  4.15  (upper).  The  graphs  show 
that  the  loss  per  bounce  is  zero  when  the  grazing  angle  is  less  than  the 
critical  angle.  Mode  1  has  the  lowest  grazing  angle  and  it  can  be  seen 
that  the  critical  angle  must  be  very  low  if  the  low  modes  are  to  suffer 
any   attenuation   at    all.      Higher    modes   suffer   higher   losses.     To  show   the 
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TABLE  4.5.   PROGRAMS  'EX.\CT'  -  REAL  AND  IMAGINARY  COMPONENTS  OF  ^^ 
AND  NUMBER  OF  ITERATIONS  TO  MEET  REQUISITE  ACCURACY. 


MODE  NUI^IBER 
2 


ACC 


REAL 
I  MAG 
ITERA.TIOMS 


REAL 
IM^.G 

ITERATIONS 


REAL 
I  MAG 
ITERATIONS 


REAL 

I  MAG 
ITEP^ATIONS 


0.212519578303 
0.000446662283 

IS 


0.  195507093057 
0.003339285283 
6 


0. 150034762533 
0. 010169974314 
5 


0.C028674''1359 
0  .037779635753 
3 


E-7   0.212512933248 
0.000450155917 
3  1 


0. 195507946339 
0.003339829695 


0.  150035442323 
0.010172092319 


0.002S57437052 
0.037779642325 


:-8   0.212512891707 
0.000450175339 
38 


0.  195507923776 
0.003339803692 
9 


0. 150035443413 
0.010172073562 
8 


0.002857437052 
0.037779642325 
d. 


E-g   0.212512888581 

0.000450175559 

CO. 


0.  195507928795 
0.003339804519 
11 


0. 150035441533 
0.010172073282 

Q 


0.002867437109 
0..037779542315 


E-10  0.212512888032 
0.000450175927 
51 


0. 195507928955 
0.003339304648 
12 


0. 150035441587 
0.010172073444 
10 


0.002357437109 
0.037779642316 
5 


E-12  0.212512387979 
0.000450176947 
58 


0. 195507923962 
0.003339804619 

14 


0.  150035441601 
0.010172073451 
11 


0.002367437109 
0. 037779542316 
5 


E-13  0.212512887974 
0.000450176949 
70 


0. 195507928961 
0.003339804620 
16 


0. 150035441602 
0.010172073449 
13 


0.002867437109 
0.037779642315 
6 
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TABLE  4.6.   PROGRAM   'EXACT'  -  REAL  AND  IMAGINARY  COMPONENTS  OF  F^ 
FOR  IMPEDANCE  BOTTOM  AS  A  FUNCTION  OF  CRITICAL  ANGLE ' 
(0  TO  .25  RADIANS) . 


MODE  NUMBER 


CRIT 
ANGLE 

1 

2 

3 

4 

5 

6 

7 

RIC- 

0 

.214751 

0 

.210105 

0.20C490 

0 

.  185134 

0 

.162416 

0 

.  128489 

0 

.063307 

BOT 

0 

.000000 

0 

.000000 

0.000000 

0 

.000000 

0 

.000000 

0 

.000000 

0 

. OCOOCO 

.25 

0 

.209522 

-0 

.209622 

0.20C611 

0 

. 185266 

0 

.162595 

0 

.  128830 

0 

.070372 

-c 

.000004 

-0 

.000004 

-0.001557 

-0, 

.002882 

-0 

.004518 

-0 

.007198 

-0 

.015854 

.20 

0 

.210226 

0 

.210226 

0.200600 

0, 

.185253 

0, 

.162589 

0 

.123325 

0, 

.070355 

-0, 

.000471 

-0 

.000471 

-0.001785 

-0. 

.003022 

.004612 

-0, 

.007258 

-0, 

.015682 

.15 

0. 

.214297 

0, 

.210212 

0.200591 

0. 

,185252 

0. 

.162584 

0. 

.  128821 

0. 

.070360 

-0. 

.000004 

-0, 

.000863 

-0.001949 

-0, 

,003126 

-0. 

,004682 

-0. 

.007303 

-0. 

.0159C2 

.10 

UNAVAILABLF 
UMA-VAILABLE 

.05 

0. 

,213628 

0. 

,213628 

0. 198850 

0. 

185216 

0. 

162580 

0. 

128792 

0. 

070C57 

-0. 

,000499 

-0. 

000499 

-0.000504 

-0. 

000591 

-0. 

002909 

-0. 

005782 

-0. 

014147 

0.0 

0. 

214649 

0. 

210107 

0.200526 

0. 

185199 
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extreme  modal  dependence  of  attenuation,  Figure  4.15  (lower)  displays 
the  loss  per  meter  for  each  mode  for  a  variety  of  cirtical  angles. 
These  curves  illustrate  that  the  increased  grazing  angles  associated 
with  the  higher  modes  result,  not  only  in  higher  loss  per  bounce,  but 
also  in  an  increased  number  of  bounces.  For  a  small  critical  angle,  say 
0.1  radians,  the  loss  per  bounce  is  approximately  three  times  as  great 
for  mode  2  as  for  mode  1,  but  the  loss  per  meter  is  about  ten  times  as 
great    (Figure    4.15). 

It  was  found  that  a  change  of  critical  angle  (change  of  the  sound 
speed  ratio)  changes  the  loss  greatly.  Figure  4.16  demonstrates  the 
change  of  pressure  due  to  mode  1  when  the  critical  angle  is  changed  from 
above  to  below  the  grazing  angle.  Changing  the  impedance  ratio  by  ten 
per  cent  caused  a  corresponding  increase  in  the  loss  but  did  not  change 
the   point   at   which  the   critical    angle   became   effective. 

To  study  the  bottom  loss  effects  for  grazing  angles  greater  than 
critical,  different  runs  were  made  for  critical  angles  that  allowed  a 
different  number  of  modes  to  propagate  in  the  low  loss  manner  (grazing 
angle  less  than  critical  angle).  Figures  4.17  and  4.18  show  how  much 
each  mode  contributes  toward  the  total  pressure,  depending  primarily 
upon  whether  its  grazing  angle  is  above  or  below  the  critical  angle. 
They  show  that  for  each  mode,  when  the  grazing  angle  is  less  than  the 
critical  angle,  the  propagation  loss  for  that  mode  is  identical  to  that 
for    the   rigid   bottom   case. 

As  a  last  demonstration  of  the  effects  of  the  impedance  bottom  on 
propagation,    data   were   compared   for   the   isovelocity    case    and    the    strong 
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Figure  4.16.   Pressure  Function  for  Mode  1  With  Critical  Angle  Above 
Grazing  (Crit  =  0.30)  and  Below  Grazing  (crit  =  0.10). 
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Propagation  Loss  for  Each  Mode  When  Cricica.I  Angle 
is  Less  Than  Grazing  Angle  for  All  Modes  (Upper) 
and  For  Mode  1  Only  (Lower). 
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and  for  Modes  1,  2,  and  3  (Lower). 


75 


positive  gradient  case  when  the  critical  angle  was  0.1  radians.  It  can 
be  seen  in  Figure  4.19  that  the  upward  refraction  casued  by  the  gradient 
lessens   loss   due   to  bottom  interactions. 


76 


RANGE  (KiM) 


Figure  4.19 


Propagation  Loss  with  Impedance  Bottom  for 
Isovelocity  Water  Conditions  (Dashed  Lines) 
and  Strong  Positive  Gradient  (Solid  Line). 
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V.      CONCLUSION 

The  program  EXACT,  which  computes  propaation  loss  by  first  solving 
for  individual  mode  wave  numbers,  was  found  to  give  realistic  values. 
The  Fast  Field  Program,  which  was  meant  to  solve  for  loss  using  an  FFT 
was  not  completed  because  of  complications  that  proved  unsolvable  in  the 
limited   time  available. 

One  of  the  problems  with  the  solution  of  the  program  EXACT,  which 
would  also  be  a  problem  with  FFP,  is  the  restriction  placed  on  the  impe- 
dance bottom.  All  of  the  sound  transmitted  into  the  bottom  is  consi- 
dered lost  to  the  bottom.  No  allowance  is  made  for  this  energy  to  be 
refracted  upward  and  re- transmitted  into  the  water.  This  might  not  be 
entirely  unrealistic  in  some  cases,  because  bottom  absorption  is  often 
very  high.  However,  this  simplification  limits  the  practicality  of  this 
program   and  is    considered    by  the   author   to   be  its   main  weakness. 

Another  point  of  consideration  is  the  number  of  modes  involved.  The 
tests  run  usually  involved  only  three  real  modes,  two  of  which  are  lost 
to  the  bottom  after  a  short  range  when  the  critical  angle  is  0.1  radians 
(sound  speed  ratio  is  1.005:1).  An  increase  of  depth  or  frequency  will 
mean  an  increased  number  of  modes.  There  might  be  U,  400,  or  40,000 
modes  present  but  a  large  portion  of  them  may  suffer  a  large  bottom 
loss.  This  also  explains  why  trapped  modes  are  so  significant  and  why 
it  is   often  sufficient   to   consider  only  a  few   modes  in  most   analyses. 


One  of  the  intended  purposes  of  the  thesis  was  to  compare  the  com- 
puter processing  time  for  each  program  to  solve  for  conditions  where  a 
large  number  of  modes  were  present.  It  was  reasoned  that  an  increase 
of  frequency  and/or  depth  would  require  a  proportionate  increase  in  the 
number  of  depth  samples  used  by  the  subroutines  that  calculated  pres- 
sure as  a  function  of  depth  for  each  horizontal  wave  number.  This  would 
mean  that  if  the  frequency  were  increased  from  50  to  100  Hertz  and  the 
depth  from  50  to  500  meters,  there  would  be  a  20  fold  increase  in  the 
processor  time  used  to  do  these  calculations  in  each  of  EXACT  and  FFP. 
However,  there  would  be  a  further  20  fold  increase  in  time  for  EXACT 
because  it  would  be  required  to  calculate  for  20  times  as  many  modes. 
The  FFP  would  not  require  any  similar  increase  and  would  therefore  have 
a  decided  advantage  when  a  large  number  of  modes  were  involved.  It  is 
possible  that  for  cases  with  a  small  number  of  modes,  EXACT  would  be 
faster,   but   the   information  was   not   available  for   comparison. 

A  simplification  in  EXACT  involves  the  accuracy  criterion  of  the  hor- 
izontal wave  number.  It  was  decided  that  errors  due  to  accuracy  limita- 
tions would  be  acceptable  if  the  wave  number  were  resolved  to  less  than 
10~^.  However,  the  subroutines  'eigens',  'modesi',  and  'modes2'  return 
the  square  of  the  wave  number.  For  the  calculated  E,^  to  remain  within 
the  limit, 

Cn    +   I  Error  I   <    ( Cn   +    10"")^ 
so 

I  Error  |  <    2  x    10"^   ^n- 
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For  the  conditions  considered,  the  correct  accuracy  criteria  should  have 
been  marginally  more  restrictive  for  all  three  modes  considered.  Had 
any  very  small  wave  numbers  resulted,  the  allowable  accuracy  criteria 
should  have  been  much  more  restrictive.  A  more  complete  program  would 
calculate  each  wave  number  to  its  own  allowable  accuracy  limit  using 
formula   5.1 

While  it  is  recognized  that  program  results  would  in  no  way  alter 
the  facts  of  nature,  the  correspondence  between  observed  phenomena  and 
predicted  results  did  add  to  the  credibility  of  the  theory  and  program 
outlined  in  this  thesis.  Many  of  the  results,  as  illustrated  in  the 
graphs,  would  prove  useful  as  teaching  aids  to  elucidate  the  principles 
of  sound  propagation  in  the  ocean.  Without  reference  to  the  development 
of  normal  mode  theory,  it  is  possible  to  illustrate  the  consequences  of 
changing  depth,  range,  frequency,  water  depth,  and  sound  speed  profiles. 
For  example,  the  program  could  be  used  to  show -the  effects  of  surface 
duct  on  long  range   propagation. 

Finally,  it  is  realized  that  the  programming  methods  used  are  far 
less  efficient  than  many  readily  available  programs.  The  intention  of 
this  thesis  was  to  develop  an  individual  Normal  Mode  program  from  a  the- 
oretical basis  without  direct  reference  to  other,  more  sophisticated  pro- 
grams. A  more  complete  study  would  compare  program  results  to  actual 
observed  data  as  well  as  to  results  of  other,  more  authoritative 
predictions . 
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APPENDIX  A 
MATRIX   METHOD   FOR   FINDING   MODAL   HORIZONTAL   WAVE   NUMBERS 

The  matrix  method  (Ref.8)  employed  in  the  program  EXACT  yields  a 
complete  set  of  eigenvalues  in  a  single  pass  utilizing  IMSL  routine 
EQRTIS.  It  is  convenient,  fast  and  has  the  advantage  that  all  modes  are 
assuredly  determined,  no  matter  how  close  together  they  are.  Accuracy 
is  not  attempted;  that  is  assured  with  the  Runge-Kutta  routines  that  are 
used  subsequently. 

The  Matrix  solution  is  a  means  of  solving  for  the  eigenvalues  of  the 
one   dimensional    differential   equation: 


^    -   ^z^'P   =  0  A.I 


where: 

((>  is  the  eigenf unction, 

k^  is   the   vertical    wave   number   i/(k^  g^) , 

^  is   the   horizontal    wave  number,  and 
z  is   the   depth. 

It   solves    the   equation   by   finding   the   eigenvalues    of    the    matrix  M    in 
the  following  Matrix  equation: 

mip   =   -(Az)2  E4)  A. 2 
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This    equation    results    from    writing    the    differential    equation    as    a    finite- 
difference   equation   at    each  of  N    depth   Az   apart.    ^^  a   column   matrix   with 
N   elements ,• each   element   defined   by: 


i+j^   =    4j(nAz)  n   =    1  ,    2 N 


A. 3 


and   M    is    a  square  symmetric   tridiagonal    matrix   which   allows    a    computation- 
ally  efficient   means   of   finding  the   eigenvalues. 


M    = 


where: 


Dx 

-1 

0 

-1 

D2 

-1 

0 

-1 

D, 

-1 


-1 


-1 


A.il 


Dn   =    2   -    (Az)^k^(nAz) 


n   =    1,    2, ,N 


A. 5 
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APPENDIX  B 
COMPUTER  PROGRAM  EXACT 


C  MAIN  PROGRAM  READS  UP  TO  30  PAIRS  OF  SOUND  VELOCITY  INFORMATICS 

C-  IT  ASSIGNS  CALCULATED  VALUES  OF  V/AVENUMBER  FOR  EACH  INCREMEKTA: 

C-  DEPTH  AND  CALCULATES  THE  NUMBER  OF  REAL  EIGENVALUES  AS  WELL  AS 

C*  MAX  WAVENUKBER  AND  THE  DEPTH  INCREMENT 

C*  ABSORPTION  IS  ASSUMED  -TO  BE  CONSTANT  AT  ALL  DEPTHS 

C'  NON-COMPLEX  EIGENVALUES  ARE  CALCULATED  BY  A  MATRIX  METHOD  IN' 

C*  SUBROUTINE  EIGENS  AND  THESE  ARE  USED  AS  A  FIRST  ESTIMATE  IN  ON'E 

C'  OF  MODES  SUBROUTINES 

C-^  PRESSURES  AND  THE  WRONSKIAN  FROM  MODES  ARE  USED  TO  CALCULATE 

C  COMPLEX  INDIVIDUAL  MODAL  PRESSURES  V7HICH  ARE  SUMMED  FOR 

C*  PARTICULAR  RANGES 

C*  BOTTOM  ABSORPTION  NOT  ACCOUNTED  FOR 

C* 


r  *  X  ■ 


COMPLEX  K( 10001) ,P( 10001) ,CALC(50) , ARGUE , VAL, CORR, INT, PDI FF , FTX , 
1PR:<,DEN0M,PPRESS(50)  ,EIG(50)  ,  PRESS  (50)  ,PRESSR,  IMAG,  ALFA,  AIG  (  50  )  , 
1FACT,REF,C0MP 

REAL  D(30) , KAY (30) , RATE (30) ,C(30),EIG2(50) , MAG (50) ,KMAX,K1, 
IKAE ( 10001 ),KAA( 10001) , PR ( 50 ) , P IM ( 50 ) 

LOGICAL  TOP 
C 

TOP  =  . FALSE . 

NUM=1 

PI  =  3. 141592654 

cma:-:=o 

CMIN=200 

FREO  =50 
DEPTHR  =  10 
DEPTHT  =37 
RMIN  =  -20 
INC  =  400 
JUMP  =  INC/ 100 
ERRMAX  =  .00003001 
CORN  =  2 
CRITA  =  . 15 
RATIO  =2.0 
12    RE.AD  (  3  ,  -  )  D  (  NUM  )  ,  C  (  NUM  ) 
C 

C*     EACH  SOUND  VELOCITY/DEPTH  PAIR  HAS  A  WAVENUMBER  CALCULATED 

c 

IF  (D(NUM)  .LT.  0)  GOTO  14 

KAY (NUM)  =  2  -  PI  -  FREQ  /  C ( NUM ) 

IF  (  C  (  NUM  )  .  GT  .  CMAX  )  THEti 

CHAX=C ( NUM ) 
END  IF 
IF  (C(NUM) . LT.CMIN)  THEN 

KMAX=KAY(NUM) 

CM IN  =  C(HUM) 

END  IF 
IF  (NUM  .Ep.  1)  GOTO  13 
C 

RATE  (NUM-1)  =  (KAY (NUM)  -  KAY(NUM-l))  /  ( D ( NUM ) -D( NUM-1 ) ) 


EILK:     ZXACT     rOi^TRAN   Al 

PRINT  103,D(r:UM-l)  ,C(NUM-1  )  ,KAY(NUM-1)  ,  RATE  (  MUM- 1  ) 
C 

13  NUM  =  NUM  +  1 
GOTO  12 

14  B0TT0M=D(NUM-1) 

N  =  IFIX( (4*B0TT0M-FRE0/CMAX  +  l)/2.0) 
V;RITE(-;,  &3S9)  N 
5339  FORMAT  (//'   THERE  ARE  ',15,'  REAL  MODES'/) 
M  =  1  +  N 

DELZ  =  BOTTOM/ (INC*N) 
PRINT  103,D(NUM-1)  ,C(I^'M-1) 

WRITE  (4, 103)D(NUM-1) ,C(NUH-1) ,KAY(NUM-1) 
IF  (C(NUM-l) . LT.C( 1) )  TOP  -     .TRUE. 
IF  (M.GT.50)  THEN 
WRITE (4, 100)  N 
100      FORMAT('TOO  MAMY  MODES  FOR  THIS  PROGRAM!!    H(REAL)  =  ',14) 
GOTO  11 
END  IF 
C 

C* 

C*  THE  NEXT  SECTION  CALCULATES  THE  EXACT  WAVENUM3ERS  FOR  THE  SU3- 

O*  ROUTINES  EI GENS,  MODES 1  AND  M0DES2 

C*  THE  REAL  COMPONENTS  FOR  MODSSl  AND  M0DES2  ARE  CALCULATED  FOR 

C*  INC*N  +  1  DEPTHS  STARTING  WITH  DEPTH  1  AT  THE  SURFACE  AND  INC'^^N-' 

C*  AT  THE  BOTTOM 

C*  AT  REGULAR  DEPTH  INCREMENTS  A  REAL  VALUE  OF  K  IS  TkKEU    FOR     -^ 

C*  USE  IN  SUBROUTINE  EIGENS       THE  NUMBER  OF  SAMPLES  IS  ICO 

C* 

C 

KOUNT  =  C 
■  NUM  =  1 
DO  18  IG  =  1, INC*N+1 
DP  =  ( IG-1)  *  DELZ 
KAA ( I G )  =  KAY ( NUM )  +  RATE ( NUM ) * ( DP -D ( I  iUM ) ) 

C      IT  IS  POSSIBLE  TO  INSERT  A  LOOP  TO  CHOOSE  ONLY  A  FRACTIOr,'  OF 
C      THE  DEPTH  SAMPLES  OF  K(  )   A  TEST  RUN  WAS  ::ADE  TO  COMPARE  THE 
C      VALUES  OBTAINED  FOR  VARIOUS  'JUMPS' 

IF  (INC.LE.200)  JUMP  =  2 
IF  (MOD  (IG-1, JUMP)  .EQ.  0)  THEN 
KOUNT  =  KOUNT  +  1 
IF(IG.EQ.l)  THEN 
KOUNT  =  0 
GOTO  886 
END  IF 

KAE( KOUNT)  =  KAA( IG) 
885      END  IF 

IF  ( DP . GE . D ( NUM  +  1 ) )  NUM  =  NUM  +  1 
18    CONTINUE 
C 

CALL  EIGENS(  BOTTOM,  KCU?IT,  M,  KAE,EIG2,  RMIN) 
C 
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KXACT  FOKIR/.N 


C-  ABSORPTION    IN    NEPERS/'METER    CALCULATED    AS    FUNCTION    OF    FREQ 

FREK    =    FREQ/1000    • 
FREEK    =    FREK'^-2 

ALFAW    =     .00001 -FREEK/ (1+FREEK) 
C 

C*  NEXT    FOLLOWS    A    LOOP    WHICH    CAN    EE    USED    TO    DETERMINE    THE    PRESSURE 

C*  OR    ATTENUATION    AT    ALL    DEPTHS    FROM    TO?    TO    BOTTOM 

C 

DO    54    IDR    =    37, 37 
DEPTHT    =    i.O-IDR 
DO    55    I    =1,M 
CRIT=CRITA 
NC    =    1 
ALFAB    =0.0 
ALFAT    =    ALFAW    +    ALFAB 
ALFA    =    CMPLX{ 0.0, ALFAT) 
5500  DO    598    IK    =    1,INC*N+1 

K( IK)    =    KAA( IK)     +    ALFA 
598  CONTINUE 

ARGUE=K(20)*-2- (PI-( I- . 5 ) /BOTTOM ) * *2 
CALC(I)    =    CSQRT( ARGUE) 

C*  CALC    IS    USED    TO    EVALUATE    THE    CORRECTNESS    OF   THE    PROGRAM 

C*  BY    INPUTTING   AN    ISOVELOCITY    PROFILE    IT    IS    POSSIBLE    TO    CALCULATE 

C*  THE    CORRECT    HORIZONTAL    WAVENUMBER    FOR    EACH    MODE    AND   THEN   COMPARE 

C*  THE    RESULTS    OF    THE    PROGR.".H    TO    THIS     'CORRECT'     ANSWER 

r^ic  ie  ic  -fr  iz  7:  ^  JT  -k  -k  i^  if  -k  iv  -k  ic  -k  -K  "K  -k  Tr  "k  yc  -k  -k  "K  ii:  ir  -i:  -k  -x  -ff  -k 

VAL    =    EIG2( I ) 
CORRN    =    1  - 
C 

IF    (TOP)    THEN 

C*  IF    THE    SOUND    VELOCITY    IS    GREATER    AT    THE    TCP    THAN    AT    THE    BOTTOM 

C*  THEN    SUBROUTINE    MODESl    IS    USED;    OTHERWISE    M0DES2    IS    USED 

C*  AFTER    EACH    RUN    OF    EITHER    SUBROUTINE,    THE    CORRECTION   MADE    TO    THE 

C*  SQUARE    OF    THE    HORIZONTAL    i.'AVELENGTH    IS    COMPARED    TO    A    ERROR 

C*  CRITERION;     IF    IT    EXCEEDS    THE    SET    LIMIT    THE    SUEROUTIN'E    IS    CALLED 

C*  AGAIN.        IF   THE    CORRECTION    IS    LESS    THAN   THE    LIMIT,     THE    PROGRAM 

C*  PROGRESSES    AND    SOLVES    FOR    THE    PRESSURE    FOR    THAT    MODE 

C* 

180  CALL    MODESl (VAL, DEL2, N, K, P, CCRR, INT, INC, RATIO, GRIT, CORRN, I 
NC    =    NC    +    1 

C  IF(CABS(VAL) .GT.CABS(K( INC-N+1 )'-2 ) )THEN 
C  GRIT    =.5*CRIT 

C  CORRI'l    =    CORN 

C  IF    (CRIT.LT. .005)    GO    TO    66 

C  VAL    =    EIG2( I ) 

C  GOTO    180 

C  END IF 


IF(CABS(CORR) .GT.ERRMAX)     GOTO    180 
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tUr.-;v.-^rj       .-.1 

IF    (CRIT.LT.CRITA) 

TH 

CRIT    =    CRITA 

GOTO    180 

END  I F 

GOTO    57 

WRITE(4, 570) I 

TLE:    EXACT 


C 
C 
C 

C 

■C55 
C 

C  IF    (VAL.EQ.O)    THEN 

ELSE 
190  CALL    M0DES2(VAL,DELZ,M,  K,?,CORR,  INT,  I'lC  ,  RATIO,  CRIT,  CO."  RI^;,  I  ) 

NC    =    NC    +    1 

IF(CABS(VAL)  .GT.CABS(K(  INC^•^>  1  )  *'^2  )  )TKEN 
CRIT    =.5*CRIT 
CORRN   =    CORN 

IF    (CRIT. LT. .005)    GO    TO    56 
VAL    =    EIG2( I ) 
GOTO    190 
E^roIF 

I  F  (  CABS  (  CORR  )  .  GT  .  ERR.MAX  )    GOTC    190 
IF    (CRIT.LT.CRITA)    THEN 
CRIT    =    CRITA 
GOTO    190 
END  IF 
GOTO    57 

56  WRITE(4, 570) I 

570  FORMAT  (  '       NO    CONVERGE^;CE  !  !  !  !  !     Flt^    THE    ERROR',  16) 

57  ENT)    IF 

EIG( I )    =    CSQRT(VAL) 

IDTl    =    IFIX(DEPTHT/(DEL2-2 ) )-2 

IDT2    =    IDTl    +    2 

PDIFF    ^    P(IDT2)     -    P(IDTl) 

DDIFF    =    DEPTHT    -     (IDT1)-DEL2 

PTX    =    P(IDTl)    +    PDIFF*DDIF~/(DELZ-2) 

IDRl    =    IFIX(DEPTHR/(DEL2-2)  ),"2 

IDR2    =    IDRl    +    2 

PDIFF    =    P(IDR2)     -    P(IDRl) 

DDIFF    =    DEPTHR    -     (IDR1)-DEL2 

PRX    =    P(IDRl)     +    PDIFF-DDIEF.'  (~ELZ-2) 

AIG(I)    =    CMPLX(REAL(EIG(  I  )  )  ,  A5S(A::-:A3(EIG(  I  )  )  )  ) 

DENOM    =    INT*CSQRT{AIG( I )-PI/2) 

PPRESS( I )=PTX*PRX/DENOM 

C         FOLLOWING    INSTRUCTIONS    ONLY    FOR    CALC    OF    R.AYLEIGH   COEFFICIENTS 

C         DONE    FOR    RANGE    OF   CRITICAL    ANGLES    .A.ND    FOR    DIFFERENT    IMPEDANCE    RATIOS 

55         CONTINUE 
C 
C 

Q^-r  ir  -k  7r  -k  it  ir  Tz  it  it  ir  ie  -k  ■k  -k  it  -k  ■x  -k  -k  it  -k  ic  •*:  it  ■k  -k  ie  -k  -re  if  -x  -x  i^  -^  i^ 

C*  THE    FOLLOWING    LOOP    CALCULATES    THE    .ATTENU.ATION    OF   THE    SIGNAL   AT 

C*  DIFFERENT    RA.NGES    BY    SUMMING    THE    COMPLEX    PRESSURE    DUE    TO    THE 

C*  INDIVIDUAL    MODES 

C*  CAN    BE    SET    UP    TO    PRODUCE    ATTENUATION   AT    REGULAR    R.ANGE    INTERVALS 

C*  OR   AT    IRREGULAR    INTERVALS    SUCH    AS    EVERY    10    METERS    OUT    TO    ICO 

C*  METERS, THEN    EVERY    100    METERS    TO    1000    AND    SO    ON. 
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U;XACT  •TOKTI^AX      Al 


C* 

IMAG=CM?LX(0, 1) 
DO    560    IRKG    =1,1 

DO    561    NRMG    =    300,3000 
R=(  10'^^■IR^;G)'■■^,'RNG 
C  R    =    50000 

FRESSR  =  0 
STRONG  =  0 
DO    562    I    =    1,N 

PRESS( I )=PPRESS( I )*CEXP( IKAG-EIG( I )-R) /SQRT(R) 
MAG(I)=    CABS(?RESS( I) ) 
FRESSR=    PRESSR    +    PRESS(I) 
PR( I )=REAL( PRESSR) 
PIM(I)    =    AIMAG(PRESSR) 

EASE    =    ATAri(AIMAG(PRESSR)/REAL(PRESSR)  ) 
STR0NG=CABS ( PRESSR ) 
C  WRITE    ('4,3434)1,     PRESS  (  I  )  ,  PFRESS  (  I  ) 

3434  F0RMAT( 16, 4E15. 7) 

C  PRINT    3  487,MAG{  I)  ,  STR0^;G 

3  487  FORMAT(2F15.5) 

FACT    =    CSQRT(K( INC*N+1)**2    -     ( EIG ( I ) /COS ( CRITA ) ) --2 ) 
REF    =    (RATIO    -    FACT )/( RATIO    +    FACT) 
FRACT    =    CABS (REE) 

COMP    =    CSQRT(K(  INC'N+ 1  ) '^^  *2 -E  IG  (  I  )  *^2  ) /' (  2^K  (  I NC* N+ 1 ) -BOTTOM  ) 
QLOS    =    ( 1 -FRACT )*CABS( COMP) 
PKAS    =   ATAN( AI MAG (REF) /REAL ( REF ) ) 
THETA    =    ARCOS(EIG(  I  )/K(  INC-'N+l)  ) 
C 
562  CONTINUE 

510  CONTINUE 

ATTEN    =    -20-ALCG10( STRONG) 

561  CONTINUE 

560      CONTINUE 

5  4         CONTINUE 
11  STOP 

END 

SUBROUTINE    MODES 1 ( VAL, DELZ , N, K, ?, CCRR, liiT, INC , RATIO , GRIT, CORRN, I ) 
COMPLEX    K( 10001) ,P( 10001) , VAL,CORR, INT, F1,A1,B1, F2 , A2 , E2 , F3,A3, 
1E3,F4, B,A, IMAG, AP,DEN0 
C 

Qtc  -K  it  -ir  -K  -K  -k  i:  -k  -k  -k  yi  -r  -k  ir  i<-  7T  -k  ir  -k  -k  -k  -k  i:  -K  if  -r-  -x  -k  -?:  yr  i:  -k  -k  -^  -k  *  -h  -r:  -i^-  -K  -^ 

C  THIS    SUBROUTINE    USES    A    FOURTH   ORDER    RUNGE    KUTTA   TECHNIQUE 

C*  TO    CALCULATE    PRESSURE    AND    ITS    DERIVATIVE       STARTING    AT    THE    SURFACE 

C*  AND    STEPPING    DOWN    A    DISTANCE    2    TIMES    THE    DEPTH    INCREMENT    UNTIL 

C*  THE    BOTTOM    IS    REACHED 

C*  CALC    AT    EACH    L    IS    IN    FACT    FOR    THE    L+1    STEP 

C-  AT   THE    BOTTOM    THE    CALCULATED    PRESSURE,    DERIVATIVE    AND    INTEGRAL 

C*  ARE    USED    TO    DETERMINE    THE    CORRECTION   TO    BE    APPLIED    TO   THE 

C*  SQUARE    OF   THE    BEST    GUESSS    OF    THE    NORMAL    MODE    HORIZONTAL    WAVE- 
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FlLi". :  EXACT     FCRTKAM   Al 

C*     ::'JM3ER     THE  CORRECTION  IS  APPLIED  AND  THE  RESULT  IS  PASSED  EAC:< 

C*     TO  THE  HAIM  PROGRAM 

C*     IN  ADDITION  THE  PRESSURE  AT  ALL  DEPTH  INTERVALS  IS  PASSED  5ACK 

C*     THIS  IS  USED  IN  THE  CALCULATION  OF  ATTENUATION 

C**^^-S!iOULD  SAVE  ONLY  THOSE  PRESSURES  REQUIRED  NOT;  ALL  i:;C-N+l 

PI  =  3 . 141592554 
INT  =0.0 
E  =  0.0 
P(l)  =  B 

DO  18  L  =  1, INC*N  -  1,2 
C 

DP  =  DELZ*(L+1) 
Fl  =  -B  -  (K(I,)**2  -  VAL) 
Al  =  A  +  DELZ*F1 
Bl  =  B  +  DELZ-A 
C 

F2  =  -Bl*(K(L+l)*-*2  -  VAL) 
A2  =  A  +  DELZ  *  F2 
B2  =  B  +  DELZ  -Al 
C 

F3  =  -B2*(K(L+1)**2  -  VAL) 
A3  =  A  +  2*DELZ-^F3 
B3  =  B  +  2-DELZ-A2 
C 

F4  =  -B3*(K(L  +  2)*'-2  -  VAL) 
B  =  B  +  2*DELZ'^(A  +  2*A1  +  2'A2  +  A3)/6 
P(L+1 )  =  B 

A  =  A  +  2-DEL2*(Fl  +  2-*F2  +  2'^F3  +  F4)/5 
INT  =  INT  +  B*^2'^-DELZ'2 
C 
18    CONTINUE 

IMAG  =  CMPLX(O.G, 1.0) 
FILL  =  REAL(K( INC*N+1) )**2 
DENG  =  CSQRT((FILL  -  VAL) /FILL) 

AP  =  CSQRT  (  F I  LL- VAL/COS  (GRIT  )  -  '  2  )  '-^  E'^  I  MAG/RAT  I  O^DENO 
CORR  =  B*(A+AP)/INT 
VAL  =  VAL  -  CCRR/CORRN 
21    RETURN 
END 

SUBROUTINE  I-:0DES2  (  VAL,  DELZ  ,  N,  K  ,  P  ,  CORR,  INT,  n.'C  ,  R.ATIC  ,  GRIT,  CORRi\',  I) 
COMPLEX  K( 10001) ,P( 10001) , VAL, CORR, INT, Fl , Al , Bl , F2 , A2  ,  E2  , 
1F3,A3,B3, F4,B,A, IMAG,DENO 
C 

C*  THIS  SUBROUTINE  USES  A  FOURTH  ORDER  RUNGE  KUTTA  TECHNIQUE 

C*  TO  CALCULATE  PRESSURE  ANT3  ITS  DERIVATIVE   STARTING  AT  THE  BGITOM   ' 

C*  AND  STEPPING  UP  A  DISTANCE  OF  2  TIMES  THE  DEPTH  INCREMENT  UNTIL    ' 

C*  THE  SURFACE  IS  REACHED 

C*  CALC  AT  EACH  L  IS  IN  FACT  FOR  THE  L-1  STEP 

C*  AT  THE  TOP,  THE  CALCULATED  PRESSURE,  DERIVATIVE  AND  INTEGRAL 

C*  ARE  USED  TO  DETERMINE  THE  CORRECTION  TO  BE  APPLIED  TO  THE 
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C*  SQUARE    OF    THE    BEST    GUESSS    OF    THE    NORMAL    MODE    HORIZONTAL    V,'AVZ- 

C*  NUMBER  THE    CORRECTION    IS    APPLIED    AND    THE    RESULT    IS    PASSED    BA^K 

C-  TO    THE    MAIN    PROGRAM 

C*  IN    ADDITIG!!    THE    PRESSURE    AT   ALL    DEPTH    INTERVAl^S    IS    PASSED    BACK, 

C-  THIS    IS    USED    IN    THE    CALCULATION    OF    ATTENUATION 

PI    =    3. 141592654 
INT    =    0 

L FLOOR    =    INC'N    +    1 
B    =    1. 

P( INC'N)    =    B 
I MAG    =    CMPLX ( 0 . 0 , 1 . 0 ) 
FILL    =    REAL(K(LFL00R) )--2 
DENO    =    CSQRT((FILL    -    VAL)/FILL) 

A    =    CSQRT(FILL-VAL/C0S(CRIT)**2)*B*IMAG/RATI0/DEN0 
A    =    0. 

IF    ( I . GT . N )    A    =    0 . 
DO    28    L    =    1, INC-N-1, 2 
C 

DP    =    DELZ*(LFL0GR-L-2) 

Fl    =    -B    *     (K(LFL00R-L+l)*-2    -    VAL) 

Al    =    A    +    DELZ-Fl 

Bl    =    B    +    DELZ*A 
C 

F2    =    -Bl'-(is(LFL0GR-L)*'^2    -    VAL) 

A2    =    A    +    DELZ*F2 

B2    =    B    +    DELZ-Al 
C 

F3    =    -B2-(K(LFLOOR-L)-*2    -    VAL) 

A3    =    A    +    2'-DELZ*F3 

B3    =    B    ^    2''DELZ'A2 
C 

F4    =    -E3'^(K(LFLOOR-L-l)--2    -    VAL) 

B    =    B    +    2-DELZ*(A    +    2-Al    +    2-A2    +    A3)/6 

?(LFLC0R-L-2)    =    B 

A    =    A    +    2-DELZ-(Fi    +    2-F2    +    2*F3    +    F4)/6 

INT    =    INT    +    B*--2-DELZ'-2 
C 
28         CONTnrJE 

CORR    =    B^-A/INT 
VAL    =    VAL    +    CORR/CORRN 
2 1         RETURN 
END 

SUBROUTINE    EIGENS ( BOTTOM, KOUNT, M, KAE , EIG2 , RM IN) 
DIMENSION    E( 10001) ,D( 10001) ,AK2( 10001) ,EIG2(50),EIC(50) 
REAL    KAE (10001) 
C  PRINT    485  5, BOTTOM, KOUNT, M,RM IN 

Q  tIt  lir  X  *  v:  *  -x  Tlr  tV  -x  IV  7f  TT  ^t  r:  -jir  i«:  *  ■^-  T\-  -^ 

C*  '■ 

C*  COMPLEX    NOT    USED    IN    THIS    SUBROUTINE    BECAUSE    IT    CALLS    IMSL  * 

C^  SUB    EQRTIS    AND    IT    IS    IMPRACTICAL    TO    MODIFY    THAT    FOR    COMPLEX 

C*  THE    REAL    EIGENVALUES    ARE    THEREFORE    CALCULATED    AND    USED    IN  * 

C*  MODESl    AND    M0DES2    WHERE    COMPLEX    EIGENVALUES    ARE    USED  ^ 
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•TILE:    liXACT  rORTRA:; 


C*  THIS    SUBRCUTIME    IS    MOT    USED    IK    THE    FFT    PROGRAM 

C* 

DZ    =    BOTTOM/ {;<0UN"T-1) 
DZ2    =   D2*D2 
C 

DO    100    I  =  l,KOU:iT 

AK2( I)    =    KAE( I)*-2       . 

D(I)    =    -AK2(I)'D22    +2.0 

E( I )    =    1.0 
100    CONTINUE 
C 

D ( KOUNT )    =    D ( KOUMT )     -    1  .  0 

E(l)    =    0.0 

IS',-/    =    0 

lER    =    0 

CALL    EQRTIS    ( D , E, KOUNT, M , I 3W, lER ) 
C 

DO    6000    J    =    1,M 
EIG2(J)    =    -D(J)/DZ2 

C  ALTHOUGH    THE    FORMULA    IS    EXACT    THE    SUBROUTINE    RUNS    IN    SHiGLE 

C  PRECISION    SO    ONE    MUST    BE    CAREFUL 

Q  :^  X  r;  X  X  ■«■  Tf  >■  TP  •?:  :r  w  w  *■  :ir  r,'  »r  i-  •>:  »r  7:  ^ 

EIG2EX    =    Ak2(20)    -     (  3  .  1415926535- (  2- J- 1 )/(  2  .  0*BOTTOM  )  )  - '•2 
6000    CONTINUE 
C 
9000    FOR.MAT    (  I  5,  3  F.  5.  10) 
8000    RETURN 

ZV.D 
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APPENDIX   C 
COMPUTER   PROGRAM    FFP 


FILE:     E-TP  FORTR.^.M      Al 


COMPLEX    WKON,U,DU, V,DV, V1,DV1, K( 25000) , FN, E IGVAL , VAL , IMAG, FUNG, 
lAREA ,  S IMPLE  ,  SOLN ,  VEE  ,  YEW ,  VRTV/N  ,  WROMK 

REAL    KMAX,D(30) ,KAY(30) , RATE (30) ,C(3G) 
C 

MUM=1 

PI    =    3 . 14159265 

CMAX=0 

CMIN=2000 
C 

FREQ    =50 

DEPTHT    =    10 

DEPTHR    =    37 

RMIN    =    40 

C  INC    ,    THE    NUMBER    OF    INCREMENTS    PER    REAL    MODE,     MUST    BE    EVEN) 

r^ir  -k  -k  -k  -^  -K  -k  -k  ir  -k  -k  i^  -x  -fr  -k  yi:  -k  -x  -k  -h  -k  it  -k  -rr  yc  -x  -k  -k  it  -K  -x  Jr  -k  -k  ^ 

INC    =2  00 
12         READ    (3,-)    D(iIUM)  ,C(NuM) 


C 


IF    ( D ( NUM )     . LT .     0 )    GOTO    1 4 
KAY{NUM)    =    2    *    PI    *    FREQ    /    G ( NUM ) 
IF    ( C ( NUM ) . GT . GMAX )    THEN 

CMAX=C ( NUM ) 

DMAX=D ( NUM ) 
END    IF 
IF    (C(NUM) .LT.CMIN)    THEN     ■ 

KMAX=KAY(NUM) 

DKMAX=D ( NUM ) 

CM  IN    =    C(r;UH) 

END    IF 
IF    (NUM    .EO.     1)    GOTO    13 

RATE     (NUM-1)     =     (KAY (NUM)     -    KAY(NUM-l))     /     ( D ( NUM ) -D ( NUM- 1 ) ) 


13  NUM    =    NUM    +    1 
GOTO    12 

14  BOTTCM=D( NUM-1) 

N    =    I  F I X ( 2  - BOTTOM" FREQ/CMAX  + . 5  )  ■ 

INCTOT    =■ INC*N 

IF    ( INCTOT. GT. 18000)     INCTOT    =    18000 

DELZ    =    BOTTOM/ INCTOT 

IF    (N.GT.50)    THEN 

WRITE (4, 100)    N 
100  F0RMAT('T00    l-lAlVf    MODES    FOR    THIS    PROGRAM!!  N(REAL)     =     ',14) 

GOTO    11 
END    IF 
C 

Q-k  -fr  -^  -k  -k  y^-  *  -K  *  -*:  k  -k  -k  *  *  i:  -k  -k  -k  -k  it  -k  -k  i^  -^^  -k  it  -k  -K  it  i:  -k  -k  -k  -x  -k  -x  -K  i<:  -k  *  -k  iv  -f.-  -K  -K  -k  -K  y^ 

c 

C      VALUE    OF    ABSORPTION    DEPENDEiN'T    OU    FREQUENCY    ONLY-NO    TEMP    DEPENDENCE) 
C 

Qyt  k  *  -k  k  k  k  k  k  *  *  k  k  ir  *  ir  k  -k  k  k  k  *  k  k  *  k  k  *  ■!:  k  k  k  k  k  it  i.  k  k  k  *  *  k  *  *  -k  *  k  k  k  it  k  i^  *  k  *  k  it  k  it  -k  ir  -k  i:  *  k  k  k  k  k  k  -. 

FFREQ    =    FREQ/ICOO 

ALFAl    =     .0000092/(.7    +    FFREQ**2) 

ALFA2    =     .0046/(6000    +    FFREQ-'--2) 
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FILE:     :•■;."?  FORTRAN'      Ai 


ALFA3    =     .000000045 

ALFA    =    (ALFAl    +    ALFA2    +    ALFA3 )     *    FFREQ-*2 
ALFA    =    .00001 
C  ALFA    =    0 

C**x*       EXAGGERATE    ALFA    TO    SEE    IF    NUM3ER    OF    ITERATIONS    IMPORTANT 
NUM    =    1 

DO    18    IG    =    1,  r:;C*N    +    1 
DP    =    (IG-1)     *    DELZ 

K(IG)    =    KAY (NUM)     +    RATE ( NUH )*( DP-D ( NUM ) ) 
IF    (DP.GE.D(NUM    +1))    NUM    =    NUM    +    1 
18         CONTINUE 
C 

IF    (DEPTHT.GT.DEPTHR)    THEN 
DUP    =    DEPTHR 
DLVl    =    DEPTHT 
ELSE 

DLW    =    DEPTHR 
DUP    =    DEPTKT 
END  IF 
C 

lUl    =    IFIX(DUP/(DEL2'2) )*2 
IU2    =    lUl    +    2 

I  LI    =    IFIX(DLV//(DELZ-^2)  )*2 
IL2    =    ILl    +    2 
C 
C 

C  ALFA    IS    CALCULATED    FOR    ONE    FREQUENCY   AND    IS    PART    OF    VAL(COM?LEX) 

C 

c 

AREA    =0.0 
IFFP    =    2047 
DO    55    I    =1, IFFP 
C 

C 

C  DURING    THIS    PHASE    OF    THE    PROGRAM   M0DES2    PROVIDES    VALUES    FCR 

C  PRESSURE    AND    DERIVATIVE    AT    BOTH    THE    TRANSMITTER    AND    THE    RECEIVER 

C  STARTING   WITH   THE    INITIAL    CONDITIONS    AT    THE    SURFACE 

C  (IE    PRESSURE    =    0         AND   DERIV    IS    AN    ARBITRARY    VALUE) 

C     .         MODES    2    PROVIDES    VALUES    FOR    PRESSURE    AND    DERIVATIVE    FOR    THE 

C  LOWER    POINT   ONLY    STARTING   WITH   THE    INITIAL    CONDITIONS    AT    THE 

C  BOTTOM       (IE    DERIVATIVE    =    0    AND    PRESSURE    IS    AN    ARBITRARY   VALUE 

C  A   VALUE    FN    REPRESENTS    THE    CONTRIBUTION    OF    EACH    'INCREMENT' 

C  AND    THE    AREA    UNDER    THE    FN   CURVE    REPRESENTS    THE    TOTAL    'PRESSURE 

C 

Q-k: 

c 


T-k'k-k-A-ivir-k-k*^ir*ifiK'*ririr-kic-k*-k-K-x-k7r-k-k-k-k-*:yr-k-! 


R    =    500 

WVNMCH    =    1 . 1*KMAX/IFFP 

V/VNM    =    I*  WVNMCH 

EIGVAL    =    CMPLX(WViMM,ALFA) 

VAL    =    EIGVAL**2 

CALL    MODESl (VAL, DELZ, INCTOT, K, U, DU, lUl, IU2 , DUP , DLW, VI , DVl , 
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FCR 


1    IL1,IL2) 

CALL  M0DES2  (VAL,  DEL2,  INCTOT  ,  K,  V,  DV  ,  ILl ,  IL2  ,  DLVI ,  INC) 
WRON  =  V*DV1-V1-DV 
878      FORMAT  (10E8.1) 

C         CALCULATE  PARTIAL  FIELD  IMTEGKAL  FOR  PARTICULAR  VALUE  OF 
C         WAVENUMBER   COMPARE  TO  THAT  CALCED  BY  PROGRAM 

C  VRTV.'N      =    CSQRT(K(20)**2    -    VAL) 

C  YEW    =    SIN(VRTWIs*DU?) 

C  VEE    =    COS(VRTWN'^(B0TTOM-DLV;)  ) 

C         WRONK  =  VRTWM-COS(VRTWN-BOTTOM) 

C         SOLN  =  YEW  *  VEE  /  WRONK 

C         SIMPLE  =  U*V/WRON 

FN  =  U*V*CSORT{EIGVAL)/WRON 

I MAG  =  CMPLX(0, 1) 

FUNC  =  FN-^CEXPl  IMAG*EIGVAL-R)^SQRT(2/(PI-R)  ) 

CONT  =  CABS (FUNC) 

AREA  =  AREA  +  FUNC*WVNMCH 

TOTAL  =  CABS (AREA) 

IF  (MOD( I, 5) .EQ.O)  WRITE  (4,190)  E I GVAL, CONT , TOTAL 

IF(MOD(  I  ,  100)  .EQ.O)  PRirJT  190,  E  I  GVAL  ,,  CONT,  TOTAL 
55    CONTINUE 
C 
11    STOP 
END 

C  MODESl  PROVIDES  PRESS  AND  DERIV  AT  BOTH  UPPER  AND  LOWER 

C  INPUTS 

C  VAL  -  THE  It:CREMENTAL  HORIZONTAL  WAVENUM  -  UP  TO  C192  VALUES 

C  DELZ-  DEPTH  INCREMENT  FN  OF  DEPTH,  N'UMB  OF  MODES,  AND  NUMBER 

C  OF  INCREMENTS  /  MODE 

C  N  -  NUMBER  OF  MODES 

C  K  -  HORIZONTAL  '..'AVErrUMBER  (COMPLEX)  AT  EACH  DEPTH  INCKEML^NT 

C  DUP  -  DEPTH  OF  UPPER  OF  TXER  AND  RXER 

C  DLW  -  DEPTH  OF  LOWER  OF  TXER  RXER 

C  IU1/IU2  INCREMENT  NUMBER  FOR  INC  ABOVE  AND  BELOW  'UPPER' 

C  IL1/IL2  INCREMENT  NUMBER  FOR  INC  ABOVE  AND  BELOW  'LOWER' 

C  INC  -  NUMBER  OF  INCREMENTS  OF  DEPTH 

C  OUTPUTS 

C  Ul  -  'PRESSURE'  AT  UPPER 

C  DUl  -  'DERIVATIVE  AT  UPPER 

C  VI  -  'PRESSURE'  AT  LOWER 

C  DVl  -  'DERIVATIVE'  AT  LOWER 

C 

C 

SUBROUTI'IE  MODESl (VAL, DELZ, INCTOT , K , U , DU, lUl, IU2 , DU? , DLW, VI , DV: 
IILI, IL2) 

COMPLEX  K( 25000) , P ( 2 ) , DP ( 2 ) , DPI ( 2 ) , PI ( 2 ) , Fl , Al , Bl , F2 , A2 , 82 , 
1F3,A3,B3, F4,A,B,PDIFF,P1DIFF,QDIFF,Q1DIFF,U, V1,DU,DV1,VAL 
C 

B  ^    0.0 

A  =  .00001 
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tTP     fcrtr;.:;  ai 

DO  18  L  =  1, IMCTOT  +  1,2 

Fl  =  -E  *  (;:(L)*-2  -  VAL) 
Al  =  A  +  DEL2:--F1 
Bl  =  B  +  DEL2-A 

F2  =  -Bl*(:<(L+l)-*2  -  VAL) 
A2  =  A  +  DEL2  *  F2 
32  =  B  +  DEL2  *A1 

F3  =  -32' (K(L+l)--2  -  VAL) 
A3  =  A  +  2-DELZ-F3 
B3  =  B  +  2'DEL2*A2 

F4  =  -B3*(K(L+2 )*-2  -  VAL) 

B  =  B  +  2-DELZ*(A  +  2*A1  +  2*A2  +  A3)/6 

A  =  A  +  2-DELZ-(Fl  +  2*F2  +  Z-FS  +  F4)/6 

IF( (L+1) .EQ. lUl)  THEN 

P(l)  =  B 

DP(1)  =  A 
END  IF 
IF  ( (L*l) .EQ. IU2)  THEN 

?(2)  =  £ 

DP(2)  =  A 
END  IF 
IF  ( (L+1) .EQ. ILl )  THEN 

Pl( 1)  =  B 

DPl(l)  =  A 
END  IF 
IF  { (L+1) .EO. IL2)  THEN 

Pl(2 )  =  B 

D?l(2 )  =  A 
END  IF 

18    CONTINUE 

PDIFF  =  P(2)  -  ?(  1) 

PIDIFF  =  Pl(2 )  -  Pl(l) 

QDIFF  =  D?(2)  -  DP(  1) 

QIDIFF  =  D?l(2)  -  DP1( 1) 

DDIFl  =  DUP  -  (IU1)*DELZ 

DID!  Ft  =  DLV,"  -  (ILl)-DELZ 

U  =  P(l)  +  PDIFF'DDIFF/(DSLZ*2) 

VI  =  Pl(l)  +  P1DIFF-D1DIFF/(DELZ^2) 

BU  =  DP(1)  +  CDIFF*DDIFF/(DELZ'2) 

DVl  =  DPl(l)  +  QlDIFF-DlDIFF/(DELZ-2 ) 
21    RETURN 

Er-!D 

C  M0DES2  PROVIDES  PRESS  AND  DER  AT  'LOV.'IR'  ONLY 

C  INPUTS 

C       VAL  -  THE  INCREMEriTAL  HORIZONTAL  V/AVENUM  UP  TO  8192  VALUES 

C       DELZ-  DEPTH  INCREMEtiT  FN  OF  DEPTH,  KUMB  OF  MODES,  AND  NUMBER 


C 
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"ILE:  Dl-P      fortk;'..\' 


C  OF  INCREMENTS  /  "ODE 

C  M  -  NUMBER  OF  MODES 

C  K  -  H0RI2CNTAL  WAVE^!UMEER  (COMPLEX)  AT  EACH  DEPTH  INCREMENT 

C  DUP  -  DEPTH  OF  UPPER  OF  TXER  AND  RMER 

C  DEW  -  DEPTH  OF  LCV.'ER  OF  TXER  RXER 

C  rJl/IU2  INCREriENT  N'UMBER  FOR  INC  ABOVE  AND  BELOW  'UPPER' 

C  IL1/IL2  INCREME::T  number  FOR  INC  ABOVE  AN'D  BELOW  'LOWER' 

C  INC  -  ^]UMBER  OF  INCREMENTS  OF  DEPTH 

C  OUTPUTS 

C  Ul  -  'PRESSURE'  AT  UPPER 

C  DUl  -  'DERIVATIVE  AT  UPPER 

C  VI  -  'PRESSURE'  AT  LOWER 

C  DVl  -  'DERIVATIVE'  AT  LOWER 

C 

C 

SUBROUTINE  M0DES2 ( VAL , DELZ ,  INCTOT, K, V, DV ,  I  LI,  IL2,DLW) 
COMPLEX  K( 25000) , P ( 2 ) , DP ( 2 ) , Fl , Al , Bl , F2 , A2 , B2 , 
1 F3 , A3 , B3 , F4 , A , B , PD I FF , QD I FF , V , DV , VAL 
B  =  .00001 
A  =  0.0 

DO  28  L  =  1, INCTOT- 1,2 
C 

LI  =  INCTOT  -  L  +  1 
Fl  =  -B  "     (K(Ll)''-2  -  VAL) 
Al  =  A  +  DELZ*F1 
Bl  =  3  +  DEL2-A 
C 

F2  =  -Bl*(K(Ll-l)--2  -  VAL) 
A2  =  A  +  DELE*F2 
B2  =  B  +  DELZ-Al 
C 

F3  =  -B2*(K(L1-1 )**2  -  VAL) 
A3  =  A  +  2*DELZ*F3 
B3  =  E  +  2'DELZ^A2 
C 

F4  =  -B3^(K(Ll-2)--2  -  VAL) 
B  =  B  +  2'-DELZ-(A  ^    2-Al  +  2'A2  +  A3), '6 
A  =  A  +  2^DELZ-^(Fi  +  2^--F2  +  2'F3  -^  F-l-j/S 
C 

'lF(  (Ll-2)  .EQ.  I  LI)  THEN 
P{1)  =  B 
DP(1)=  -A 
END  IF 

IF( (Ll-2) .EQ. IL2)  THEN 
P(2)  =  B 
DP(2)  =  -A 
EMDIF 

28    CONTINUE 

PDIFF  =  P(2 )  -  P(l) 

QDIFF  =  DP(2)  -  DP( 1) 

DDIFF  =  DLW  -  (IL1)^DEL2 

V  =  P(l)  +  PDIFF*DDIFF/(DELZ-2) 

DV  =  DP(1)  +  QDlFF-'DDlFF/iDELZ*!) 

21         RETURN 
END 


C 


C 
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